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Abstract 

It  is  the  purpose  of  the  proposed  research  to  develop  a  digital  computer  simulation 
model  to  study  the  effects  of  an  internal  fault  in  a  large  permanent  magnet  ac  synchronous 
motor.  Permanent  magnet  motors  are  being  considered  as  an  alternative  for  ships  with 
electric  propulsion  systems.  In  an  electric  propulsion  system  a  large  motor  will  be  directly 
connected  to  a  propulsion  shaft.  A  windmilling  shaft  will  continue  to  turn  the  rotor  of  the 
propulsion  motor  after  the  motor  has  been  disconnected  from  its  electrical  power  supply 
source. 

Following  an  internal  electrical  fault  in  a  propulsion  motor,  it  is  expected  that  the 
motor  will  be  disconnected  from  its  electrical  supply  source.  With  the  ship  operating  at  or 
near  rated  speed  following  a  casualty  to  the  propulsion  plant,  the  ship  will  coast  down  to  a 
stop  or  until  the  crew  takes  action  to  stop  the  ship.  A  windmilling  permanent  magnet 
motor  will  generate  a  large  enough  internal  voltage  to  continue  to  support  large  fault 
currents. 

This  research  will  focus  on  the  fault  transient  and  the  motor  behavior  during  the 
time  that  the  propulsion  shaft  is  windmilling.  Shorting  the  motor  terminals  will  be 
considered  as  a  means  of  reducing  the  power  input  into  the  fault. 
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Chapter  1.  Introduction 

The  use  of  electric  ship  propulsion  offers  significant  advantages  in  ship  design, 
construction  and  operation.  Placing  electric  propulsion  motors  as  far  back  in  the  ship  as 
possible  serves  to  reduce  long  shaft  lengths.  Propulsion  shafts  often  go  through  many 
compartments,  creating  design  complications  with  shaft  line  component  alignment  and 
compartment  arrangement.  In  addition,  electric  propulsion  can  reduce  the  number  and 
type  of  prime  movers  in  the  ship.  Propulsion  power  and  shipboard  electrical  power  can  be 
derived  from  common  prime  movers. 

In  electric  propulsion  ships,  prime  movers  can  be  located  anywhere  in  the  ship. 
This  flexibility  can  improve  survivability  and  make  maintenance  and  shipboard 
arrangements  easier  [1]. 

Using  electric  drive  can  provide  the  ship  with  greater  operational  flexibility. 
Electric  drive  ships  can  operate  the  prime  movers  at  their  optimal  speed  and  most  efficient 
speed.  By  operating  the  prime  movers  at  their  most  efficient  speed,  the  ship's  fuel 
consumption  can  be  lowered.  This  could  decrease  the  frequency  of  refueling  and  increase 
the  endurance  of  the  ship. 

Electric  power  propulsion  has  been  used  in  ships,  commercial  and  military,  for 
over  50  years.  However,  some  drawbacks  have  been  higher  costs  and  greater  weight  than 
mechanical  drive  systems.  Advanced  systems  using  permanent  magnets  have  the  potential 
of  being  lighter  than  similar  conventional  systems. 


The  US  Navy  is  developing  modern  electric  drive  propulsion  systems  for  its  ships. 
In  some  of  the  conceptual  designs  for  these  electric  propulsion  systems,  the  use  of 
permanent  magnet  propulsion  motors  and  ship  service  generators  has  been  considered. 
An  issue  that  needs  to  be  investigated  is  the  behavior  of  such  machines  during  and 
following  an  electrical  fault  inside  the  motor  or  generator,  where  the  machine  supply 
breakers  might  not  be  able  to  prevent  damage.  Of  concern  are  internal  arcing  faults  that 
can  result  in  very  large  currents  and  high  localized  temperatures  [2].  These  large  currents 
and  high  temperatures  can  cause  extensive  damage  to  insulation  and  conductors  and  more 
importantly  to  the  permanent  magnet  themselves. 

The  arcing  process  itself  can  cause  significant  damage  to  the  machine.  Arcs  can 
compromise  the  electrical,  as  well  as,  the  mechanical  integrity  of  the  machine.  Localized 
damage  caused  by  the  electric  arc,  such  as  pits  or  hardened  spots,  can  serve  as  crack 
initiators  from  which  cracks  that  can  be  induced  by  machine  vibrations  can  initiate  and 
propagate  [3]. 

After  a  permanent  magnet  generator  or  motor  is  disconnected  from  the  rest  of  the 
system  it  will  continue  to  generate  an  internal  voltage,  Eaf,  until  the  field  comes  to  a  stop. 
In  wound  field  machines  the  field  circuit  breaker  can  be  tripped  simultaneously  with  the 
main  circuit  breakers  thus  essentially  eliminating  the  field.  Figure  1  shows  a  simple 
schematic  of  an  ac  synchronous  machine. 
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Figure  1-1  Simplified  Motor 
Model 


The  internal  voltage  generated  in  the  machine  is  proportional  to  the  rotational 
speed  of  the  field.  In  big  machines,  such  as  those  used  for  propulsion  or  power 
generation,  these  internal  voltages  can  be  very  large.  The  generated  internal  voltage  can 
be  large  enough  to  continue  the  arcing  process  and  cause  further  damage  while  the 
machine  is  coasting  down. 

The  purpose  of  this  research  is  to  develop  some  tools  that  can  be  used  to 
investigate  what  happens  to  a  large  permanent  magnet  motor  following  an  internal  fault 
while  the  motor  coasts  down. 

1.1   Ship  Propulsion  Systems 

Current  naval  propulsion  systems  consist  of  multiple  diesel  engines  or  gas  turbines 
coupled  to  one  or  more  propulsion  shafts  through  a  set  of  reduction  gears  and  clutches. 
Using  two  engines  per  shaft  provides  redundancy  and  continuity  of  propulsion.  In  the 
case  of  damage  or  maintenance  to  one  engine  while  the  ship  is  underway,  propulsion  can 
still  be  maintained  on  that  shaft. 


Although  mechanical  systems  are  very  simple  and  reliable,  they  have  some  major 
disadvantages.  Mechanical  drive  systems  require  separate  prime  movers  for  electrical 
power  generation.  Alignment  between  the  shaft  and  the  prime  mover  requires  that  prime 
movers  be  located  as  low  in  the  ship  as  possible.  This  requirement  results  in  large 
amounts  of  "lost  volume"  inside  the  ship  and  the  superstructure  for  intake  and  exhaust 
trunks.  Some  of  this  volume  can  be  recovered  in  ships  with  electric  drive. 

Relatively  light  prime  movers,  such  as  gas  turbines,  can  be  located  higher  in  the 
ship  as  electric  drive  does  not  require  alignment  between  the  propulsion  motor  and  the 
prime  mover.  Conceivably,  prime  movers  could  be  located  in  the  superstructure. 
Locating  the  prime  movers  higher  in  the  ship  can  minimize  the  arrangeable  volume  lost  to 
long  exhaust  and  intake  trunks.  Other  advantages  and  drawbacks  of  these  systems  are 
discussed  in  [4]. 

Another  possibility  that  electric  propulsion  offers  is  the  capability  to  move  the 
propulsion  motors  outside  the  hull  of  the  ship. 

1.2  Permanent  Magnet  Motors 

The  machine  used  in  this  research  is  a  permanent  magnet  ac  synchronous  motor. 
This  is  essentially  an  ac  synchronous  motor  in  which  the  field  windings  have  been  replaced 
by  permanent  magnets.  The  analysis  used  for  this  machine  is  that  used  for  wound  field  ac 
synchronous  machines  assuming  that  the  machine  is  excited  by  a  field  current  of  constant 
value. 

Permanent  magnet  machines  have  a  number  of  advantages  over  wound  field 
machines.  One  of  their  most  significant  advantages  is  the  elimination  of  the  field  windings 


and  the  power  losses  generated  in  these  windings.  In  addition,  since  no  electrical 
connections  are  required  to  supply  the  field,  this  eliminates  the  need  for  slip  rings  and 
brushes.  Therefore,  eliminating  the  field  windings  can  also  result  in  smaller  machines  than 
wound  field  machines.  This  is  especially  important  in  marine  applications  that  are  volume 
limited  where  space  is  critical. 

Some  of  the  limitations  of  permanent  magnet  machines  come  from  the  permanent 
magnets  themselves.  Excessive  currents  in  the  motor  windings  or  excessive  heat  could 
result  in  demagnetization  of  the  magnets.  Other  limitations  and  advantages  of  permanent 
magnet  machines  are  discussed  in  detail  in  [5]  and  [6]. 

The  control  aspects  and  the  electrical  power  requirements  of  permanent  magnet  ac 
machines  are  not  addressed  in  this  research.  Such  aspects  of  permanent  magnet  and  other 
electrical  machines  are  addressed  in  references  such  as  [7]  and  [8]. 

For  the  purposes  of  this  research  the  permanent  magnet  motor  will  be  assumed  to 
have  been  operating  in  steady  state  at  the  time  the  fault  is  initiated.  After  the  fault  is 
initiated  the  motor  will  be  disconnected  from  the  rest  of  the  network  and  allowed  to  coast 
down.  This  research  will  only  address  large,  low  speed  motors  such  as  those  that  would 
be  used  for  ship  propulsion.  In  this  configuration  the  propulsion  shaft  will  be  directly 
connected  to  the  rotor  of  the  machine. 

The  permanent  magnet  motor  used  in  this  research  is  a  notional  machine.  A 
notional  machine  was  used  as  this  type  of  propulsion  machinery  is  not  yet  in  use  in  naval 
ships.  Design  and  sizing  calculations  of  this  machine  are  discussed  in  chapter  two. 
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One  of  the  advantages  of  electric  drive  propulsion  is  the  ability  to  directly  drive  a 
propulsion  shaft  without  the  need  for  reduction  gears.  It  is  desirable  that  ship  propulsion 
motors  operate  at  low  speeds,  so  these  machines  will  have  a  large  number  of  poles.  As 
the  number  of  poles  increases  the  physical  size  of  the  machine  will  increase.  During  the 
notional  design  phase  an  effort  was  made  to  keep  the  size  and  weight  of  the  motor 
comparable  with  current  propulsion  machinery.  For  this  research  a  flux  concentrating 
permanent  magnet  machine  was  selected.  The  geometry  of  this  machine  is  shown  and 
discussed  in  more  detail  in  chapter  two. 

1.3  Ship  Model 

Following  the  fault  to  the  motor,  after  the  main  supply  breakers  to  the  motor  are 
open,  the  ship  is  assumed  to  coast  down  to  some  fraction  of  its  initial  speed.  During  this 
coast  down  period  and  during  the  initial  phases  of  the  casualty  it  will  be  assumed  that  no 
action  is  taken  by  the  ship's  crew  to  slow  down  or  stop  the  shaft.  The  basis  for  this 
assumption  is  that  during  the  initial  phase  of  the  casualty,  the  first  indication  available  to 
the  operators  will  be  the  tripping  of  the  main  supply  breakers.  Once  the  main  supply 
breakers  trip,  there  will  exist  an  inherent  time  delay  before  action  is  taken  to  stop  the 
affected  shaft  and  take  corrective  actions.  This  delay  is  due  to  the  time  that  it  will  take  for 
ship's  personnel  to  evaluate,  recognize  and  act  to  combat  the  casualty. 

Automatic  protection  systems,  that  operate  when  the  breakers  trip,  could  be  used 
to  minimize  the  time  delay  in  stopping  the  shaft.  A  problem  with  these  systems  is  that  they 
could  initiate  protective  action  during  spurious  breaker  trips.  During  these  trips  protection 
is  not  necessary  and  could  be  detrimental  to  the  operation  of  the  ship. 

11 


The  mechanical  energy  supplied  to  the  motor  while  the  ship  is  coasting  down  is 
proportional  to  the  speed  of  the  shaft  squared.  It  can  be  shown  that  for  a  given  propeller 
the  rotational  speed  of  a  windmilling  shaft  is  proportional  to  the  speed  of  the  ship.  As  the 
ship  coasts  down  the  mechanical  energy  supplied  to  the  motor  will  decrease  as  the  ship 
slows  down. 

1.4  Faults  in  Permanent  Magnet  Ship  Propulsion  Motors 

The  focus  of  this  research  is  investigating  the  performance  of  large  permanent 
magnet  motors  such  as  those  that  would  be  used  for  ship  propulsion.  It  can  be  expected 
that  part  of  the  electrical  protection  of  these  large  motors  will  be  main  supply  breakers. 
These  breakers  will  disconnect  the  motor  from  its  electrical  power  supply  upon  detection 
of  a  fault  in  the  motor. 

Internal  faults  are  of  special  interest  due  to  the  large  amounts  of  energy  that  can  be 
supplied  to  these  motors  by  a  free  spinning  shaft  after  the  electrical  supply  breakers  are 
tripped.  Since  the  field  of  the  machine  is  supplied  by  the  permanent  magnets,  the  motor 
with  the  free  spinning  shaft  will  behave  like  a  generator  supplying  power  to  the  fault. 
Referring  to  figure  1-1,  the  excitation  voltage,  Eaf,  is  proportional  to  the  flux  produced  by 
the  field,  (f)f,  and  the  frequency  of  the  field  excitation,  a>  [5] 

Eaf=K©(j)f  (1-1) 

In  permanent  magnet  machines  a  constant  field  flux,  <J)f,  is  supplied  by  the 

permanent  magnets.  Therefore  to  stop  this  generator  action,  the  shaft  must  be  stopped. 

Stopping  the  shaft  could  be  accomplished  by  using  a  mechanical  shaft  brake,  changing  the 
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pitch  of  the  propeller  in  ships  with  controllable  pitch  propellers  or  stopping  or  slowing 
down  the  ship.  Ships  with  more  than  one  propulsion  shaft  can  use  the  unaffected  shaft  to 
slow  down  or  stop  the  ship. 

Other  ways  to  slow  down  or  stop  permanent  magnet  machines  are  discussed  in  [9]. 
These  methods  are  for  unfaulted  machines  whose  shafts  are  not  being  driven  by  an 
external  source  such  as  a  ship's  windmilling  propeller.  The  windmilling  propeller  and  shaft 
will  rotate  at  a  speed  that  is  proportional  to  the  ship's  speed. 

1.5  Research  Approach 

This  research  studied  the  dynamic  behavior  of  permanent  magnet  ac  synchronous 
motors  following  an  internal  fault.  The  motor  studied  was  assumed  to  be  directly 
connected  to  a  ship's  propulsion  shaft.  Following  the  fault,  the  motor  was  disconnected 
from  its  electrical  supply  source  and  allowed  to  windmill  as  the  ship  coasted  down  to  a 
fraction  of  its  initial  speed. 

Different  from  wound  field  machines,  permanent  magnet  machines  have  a  constant 
field  flux  supplied  by  the  magnets.  As  the  shaft  windmills  this  constant  field  flux  will 
generate  an  internal  voltage  in  the  machine.  A  sufficiently  large  internal  voltage  can 
continue  to  support  an  internal  fault  in  the  machine  as  the  shaft  windmills. 

To  accomplish  the  goals  of  this  research  the  following  tasks  were  identified  and 
performed: 

1.  A  notional  permanent  magnet  motor  having  the  performance  requirements  of  a 
ship's  propulsion  motor  was  derived.  To  do  this  task,  a  motor  design  spreadsheet  was 


13 


developed.  By  specifying  the  principal  motor  requirements  the  motor  parameters  were 
estimated. 

2.  A  dynamic  model  of  the  permanent  magnet  motor,  which  incorporated  the 
parameters  of  the  notional  design  was  derived. 

3.  Models  for  the  internal  fault  and  of  the  ship  were  derived. 

4.  The  models  of  the  motor,  fault  and  ship  were  incorporated  into  a  dynamic 
simulation  program. 

5.  The  simulation  was  run  and  the  results  evaluated. 
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Chapter  2.  Conceptual  Design 

This  chapter  describes  the  procedure  used  to  design  the  notional  permanent 
magnet  motor  used  in  this  research.  The  motor  designed  incorporates  desired  attributes  of 
a  motor  for  naval  propulsion.  The  procedure  described  in  this  chapter  is  implemented  into 
an  Excel  [10]  spreadsheet,  Appendix  A.  For  the  motor  design,  a  set  of  requirements  was 
established  and  some  initial  assumptions  were  made  as  to  the  geometry  and  operating 
parameters  of  this  motor.  The  initial  assumptions  and  performance  requirements  are 
summarized  in  Table  1 . 

Table  1.  Motor  Specifications 


Number  of  Phases 

3 

Frequency  (Hz) 

60 

Rotor  speed  (rpm) 

200 

Rated  power  (Hp) 

40000 

Operating  voltage  (V) 

1000 

Power  factor 

0.8 

Winding  factor,  kw 

0.9 

L/D 

0.22 

Tooth  fraction,  ^p 

0.5 

A  line  frequency  of  60  Hz  was  selected  for  the  motor  design  since  this  is  the 
frequency  commonly  used  in  U.  S.  Navy  ships  electric  service  applications. 
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The  type  of  machine  that  was  used  in  this  research  is  a  flux  concentrating 
permanent  magnet  machine.  In  this  type  of  machine  the  magnets  are  oriented  so  that  their 
magnetization  is  azimuthal  as  shown  in  figure  2-1. 


Wm 


Magnet 


Figure  2-1.  Motor  cross  section 

Using  the  geometry  described  in  figure  2-1  and  the  requirements  specified  in  table 
1,  some  of  the  initial  sizing  considerations  were  started.  For  the  required  rotor  speed,  n, 
and  electrical  frequency,  f,  the  number  of  pole  pairs,  p,  in  the  machine  was  determined 
using  p  =  60f/n.  A  rotor  speed  of  200  rpm  was  selected  for  this  design  because  this  speed 
is  consistent  with  rated  shaft  speeds  for  naval  ships. 

Since  this  machine  will  be  used  for  ship  propulsion,  one  of  its  desired  attributes  is 
that  it  is  comparable  is  size  and  weight  to  conventional  propulsion  machinery.  By 
calculating  the  radius  and  the  length  of  the  machine,  an  initial  determination  of  the 
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feasibility  of  the  motor  can  be  obtained.  Weight  calculations  are  performed  once  all  the 
machine  components  are  sized. 

Once  the  rated  speed  of  the  rotor  was  established,  the  radius  of  the  rotor,  R,  was 
determined  from  the  mechanical  power  output  required  and  the  average  gap  shear 
stress,  (t).  To  calculate  the  radius  of  the  rotor  the  following  two  equations  were  used 

Power  =  (Torque)co  m  (2- 1 ) 

Torque  =  (t)(2tcRL)R  (2-2) 

The  aspect  ratio  indicated  in  table  1  was  calculated  using  the  following  relation 
suggested  by  Levi  in  [10] 

--V*       (2-3) 
D      2 

The  aspect  ratio  estimated  using  equation  (2-3)  minimizes  the  winding  resistance  for  a 

given  electromotive  force. 

Using  the  aspect  ratio,  L/D,  calculated  from  (2-3)  and  combining  equations  (2-1) 

and  (2-2),  the  radius  of  the  rotor  is 


R  =  3         P°*;7      ,  (2-4) 

Once  the  machine  radius  was  calculated,  the  active  length  of  the  machine  was  determined 
as  L  =  2R(L/D) . 

A  value  for  the  average  gap  shear  was  estimated  in  order  to  calculate  the  size  of 

the  machine.  The  average  gap  shear  stress  can  be  estimated  by  (x)  «  -^=B,K  ,  where  B\ 
is  the  peak  value  of  the  fundamental  magnetic  flux  density  and  K  is  the  root  mean  square 
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(rms)  value  of  the  effective  armature  reaction  current  density  [6].  Gap  shear  has  units  of 
pressure,  pascals,  (Pa). 

For  this  machine  design  the  value  selected  for  the  average  gap  shear  was  based  on 
machine  sizing  considerations  and  cooling  requirements.  By  maximizing  the  gap  shear,  the 
machine  radius  can  be  reduced,  equation  (2-4);  however,  high  values  of  shear  will  require 
additional  cooling  requirements. 

The  value  of  shear  stress  selected  is  that  of  an  air-cooled  machine.  This  type  of 
machine  was  selected  in  an  effort  to  maintain  simplicity.  For  this  type  of  machine,  and  for 
the  specified  rating,  a  shear  value  of  approximately  50,000  Pa  was  obtained  [6]. 

Figure  2-2  shows  the  variation  of  rotor  radius  as  a  function  of  gap  shear.  From 

equation  (2-4)  it  can  be  noted  that  R  x  —= .  Higher  shear  values  will  result  in  smaller 

VI 

machines.  However,  these  smaller  machines  will  require  additional  cooling  provisions, 
adding  to  the  complexity  and  most  likely  to  the  overall  size  of  the  machine. 


Fractional  reduction  in  Rotor  Size  as 
a  Function  of  average  gap  shear 


>/->0>/"iO>/~)0>/"10>/">Ol/~>Ol'"iO>/1 
>— i     —     <N     <N     <"*"i     fl     *3"     *3"     ^">     "~>     ^o     ^o     r-     c~- 

Average  gap  shear  (Pa  x  1000) 


Figure  2-2.  Radius  reduction  as  a  function  of  gap  shear 
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Considering  the  above  reasons  an  average  gap  shear  of  50,000  Pa  was  selected  and 
used  for  the  motor  design.  This  shear  value  is  at  the  high  end  for  typical  slow  speed  air- 
cooled  machines  [6].  For  the  specified  shear  and  aspect  ratio  the  rotor  diameter  was 
calculated  to  be  approximately  4.2  m  and  the  length  of  the  machine,  L,  approximately  one 
meter. 

2.1  Air  gap  size 

Mechanical  considerations  and  performance  requirements  drive  the  length  of  the 
air  gap,  g.  The  minimum  physical  length,  in  meters,  of  the  air  gap  can  be  estimated  from 
the  following  relation  [11] 


-3-    L 


3/4 

g*3.35xl0-^J     D  (2-5) 

Using  the  calculated  radius  of  the  rotor,  equation  (2-4),  and  L/D,  (2-3)  a  minimum 
gap  length  of  5  mm  is  obtained.  This  gap  length  was  used  for  this  research. 

2.2  Magnet  dimensions 

Considering  the  geometry  shown  in  figure  2-1,  and  assuming  that  the  magnets 
occupy  one  half  of  the  circumference  of  the  rotor,  the  width  of  a  magnet  is  given  by 

wm  =  —  and  from  figure  2-1  w     =  2Rsin— L.  Combining  these  two  equations  and 

P  2 

solving  for  the  magnet  angle,  6t 

0t=2sin'—  (2-6) 

2p 

Then  the  pole  face  angle,  0m,  is  given  by 
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e».=--et  (2-7) 


The  remaining  dimension  of  the  magnets  yet  to  be  determined  is  the  height,  hm. 
This  dimension  is  determined  by  considerations  other  than  just  geometry,  such  as  the 
magnetic  flux  density  at  the  air  gap.  In  order  to  calculate  the  magnet  height,  the  gap  flux 
density  needs  to  be  determined. 

In  flux  concentrating  machines,  the  flux  density  in  the  gap  is  greater  than  the  flux 
density  in  the  magnets.  To  calculate  the  fields  in  this  machine  a  simple  reluctance  model 
presented  in  [12]  was  used.  For  the  given  geometry,  the  flux  path  in  the  rotor  involves  a 
magnet  and  one  half  of  each  of  two  adjacent  pole  pieces.  For  one  pole  piece  the 
permeance  of  the  air  gap  is  given  by 

DO 

p„p  =  H0L— ^         (2-8) 


gap 


and  the  permeance  of  a  magnet  is 


pm  =  uL-^-  (2-9) 


m         ~o 

W. 


The  magnetic  flux  density  in  the  magnets,  Bm,  is  calculated  using  a  simplified 
magnetic  circuit.  This  circuit  consists  of  the  magnet's  own  permeance  in  series  with  one 
half  of  the  permeance  of  each  of  two  pole  pieces.  So  that  the  flux  density  of  the  magnets 
is  given  by 

B     =B   — ^EE (2-10) 

Pgap+Pm 

and  the  flux  density  in  the  gap  can  be  calculated 
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9   Y\ 

B      =B &■  (2-11) 

m  R  A 

m 

Solving  equation  (2-11)  for  the  magnet  height,  the  following  relation  is  obtained 

m  m  V  ) 

m 

For  a  given  radius  and  magnet  angle,  the  ratio  of  gap  flux  to  magnet  flux  will 
determine  the  magnet  height.  The  gap  flux  density  that  can  be  achieved  in  the  machine 
will  be  limited  by  stator  teeth  and  back  iron  saturation.  This  limit  will  be  checked  later  in 
the  design.  Once  a  magnet  height  is  selected,  the  magnet  spacing  needs  to  be  checked. 
For  the  given  geometry  the  closest  point  between  two  magnets  occurs  at  the  interior 
corners.  Using  simple  geometry  the  distance  between  two  adjacent  corners  can  be  shown 
to  be  [12] 

s=2(R-hm)sin wm  cos —       (2-13) 

V  m)       2p        m        2p 

After  checking  that  the  magnet  dimensions  are  compatible  with  the  rotor  size  and 
the  gap  flux  density  has  been  selected,  the  surface  current  density,  K,  necessary  for 

V2« 


operation  of  the  machine  at  rated  power  can  be  calculated,  K  « 


B, 


2.3  Determination  of  terminal  current  and  machine  rating 

For  a  machine  with  a  small  gap,  it  can  be  assumed  that  the  magnetic  flux  is  not  a 
function  of  radial  position.  The  space  fundamental  of  this  flux  is  then  of  the  form 

B,  =  *sin^Bgap  (2-14) 
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With  this  value  of  flux,  the  internal  voltage  can  be  estimated  using 

2RLN  km  B, 


E*  = 


s       s  1 


42 


(2-15) 


The  terminal  current  into  the  machine,  It,  is  given  by 


T      _  slots W  slots'1  slot    T 

6N. 


(2-16) 


where  Nslotswslots  =  27dU.p  and  K  =  Jahslotks.  Once  the  internal  voltage  and  the  terminal 
current  are  known,  the  rating  of  the  machine  |P  +  jQ|  =  3VtIt  can  be  determined  if  the  ratio 
of  terminal  voltage  to  internal  voltage,  v  =  Vt/Eaf  ,  is  known.  A  method  to  calculate  v 
will  be  proposed  later  in  the  chapter;  however,  this  method  requires  that  the  machine 
reactances  be  known. 


2.4  Machine  Reactances 

The  permanent  magnets  are  in  the  main  flux  path  of  the  armature.  The  presence  of 
the  magnets  will  make  the  machine  salient  since  the  direct  and  quadrature  axes  will  be 
affected  differently.  Derivation  of  the  direct  and  quadrature  axes'  reactances  is  shown  in 
several  places;  however,  for  consistency  the  notation  used  in  [12]  will  be  maintained.  The 
direct  and  quadrature  reactances,  Xa  and  Xq,  are  given  by 


Xd  = 


3  4  u0N;k;RL        .    P0r 
-V5 coy  sin — ■ 


2  71 


P  g 


3  4  u  NVRL    (       .    p0t^ 
q      2  7i        p-g  k  2  J 


(2-17) 


(2-18) 
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These  formulas  have  a  new  unknown,  the  number  of  turns  per  phase  in  the  stator, 
Ns.  The  value  of  Ns  will  be  calculated  once  v  and  the  rating  of  the  machine  are  known. 


2.5  Internal  Voltage 

To  calculate  the  ratio  of  terminal  voltage  to  internal  voltage  the  phasor  diagram 
shown  in  figure  2-3  is  used.  A  two-axis  representation  of  the  machine  is  shown  in  figure 
2-3.  The  mathematical  transformation  used  to  derive  the  two-axis  model  is  discussed  in 
chapter  three. 


d-axis 

IoK 

/       I 

/> 

v-t 

jXdld 

Iq^ 

6 

E* 

Ea 

/jXqlq 

q-axis 

Figure  2-3.  Phasor  diagram  for  Negatively  Salient 

Motor 

From  figure  2-3  the  following  set  of  phasor  relations  are  derived 


E*=^Vt2+(ltXq)2-2VtIXqsinM/ 
Id  =It  sin(\|/  +  8) 


5  =  tan 


IX  cos\|/ 

V  +  IX  cos  v|/ 


(2-19) 
(2-20) 

(2-21) 
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Eaf=E*-(Xq-Xd)ld  (2-22) 

To  solve  this  set  of  equations  an  iterative  method  is  suggested  in  [12].  The  end 
result  of  the  proposed  method  is  a  value  for  the  ratio  of  terminal  to  internal  voltage,  v.  To 
implement  the  method,  the  set  of  equations  (2-19)  through  (2-22)  and  the  machine 
reactances,  equations  (2-17)  and  (2-18)  are  normalized  with  respect  to  the  internal 
voltage,  equation  (2-13).  After  normalizing,  per  unit  scaling,  and  some  simplification  the 
reactances  are  given  by 

^RJL  V^ysin^-  (2-23) 


ad  p 


X 


Pg     '       B, 


Pg  Bj  v 


The  set  of  phasor  relations  after  normalizing  and  assuming  operation  at  rated 
current  ( it  =  1.0  p.u.)  and  power  factor  is  given  by 

e*  =  <Jv2  +xa3q  -2xaqvsinij/  (2-25) 


xaqvcos\|/ 


5  =  tan"1  — 2 —  (2-26) 

v  +  xa  vsinvj; 

id  =  sin(\|/  +  5)  (2-27) 

e«f  =  e*-(xaq  -xad)id  (2-28) 

Using  the  normalized  reactances  and  the  set  of  equations  (2-25)  through  (2-28)  the 
suggested  method  of  [12]  seeks  a  set  of  values  that  will  solve  the  above  set  of  equations  in 
which  the  fixed  point  is  eaf  =  1.0.  The  iterative  method  is  started  by  selecting  an  arbitrary 
value  for  v  and  solving  equations  (2-25)  through  (2-28)  sequentially.  For  subsequent 
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iterations  a  new  value  of  vnew  =  v^/^/eJis  used  and  the  iterative  process  is  continued 
until  eaf  converges  to  1.0.  Once  convergence  is  obtained  the  final  value  of  v  is  used  to 
determine  the  machine  rating 

3VtIt  =-^-R2IA  KB,v     (2-29) 
V2  p 

If  a  value  of  terminal  voltage  is  selected,  table  1,  then  the  machine  terminal  current 
can  be  calculated,  equation  (2-15). 

2.6  Stator  sizing 

Once  the  terminal  current  is  known  the  number  of  stator  turns  per  phase  is 
calculated,  equation  (2-16).  The  slot  width  and  height,  the  number  of  slots  and  the  size  of 
the  back  iron  need  to  be  calculated  to  complete  the  initial  sizing  of  the  machine.  If  a  tooth 
fraction  of  Xv  =  0.5  is  used  and  assuming  that  R  »g,  then  the  slot  width  is  calculated  as 

2tiR 

™s  =  -^r  (2"30) 

6N  k 


s      s 


and  the  number  of  slots  in  the  stator  is  given  by 

NsloUs=6NsksXp  (2-31) 

The  slot  height,  hs,  is  determined  principally  by  the  insulation  required  by  the 

armature  windings.  Commonly  these  requirements  are  such  that  the  copper  area  in  the 
slot  is  between  40%  to  60%  of  the  area  slot  [6].  For  this  design  the  slot  copper  fraction, 
kcm  selected  was  0.5.  With  these  assumptions  the  slot  height  is  calculated  by 

K 


s'ot  ~~  Jakskcu 
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(2-32) 


Limits  on  the  current  density,  Ja,  are  established  by  the  maximum  temperature  rise, 
6W,  in  degrees  Kelvin,  above  ambient  temperature  (40  C),  allowed  in  the  copper.  The 
allowable  temperature  rise  is  based  on  the  class  of  insulation  used  in  the  machine.  For  an 
air-cooled  machine,  the  maximum  current  density  based  on  thermal  considerations  can  be 
estimated  by 

Ja  <  h6jc"ks  (2-33) 

K. 

where  ycu  is  the  conductivity  of  copper  and  h  is  the  overall  heat  transfer  coefficient. 

Equation  (2-33)  represents  the  energy  balance  between  the  heat  generated  in  the  copper 

and  the  heat  removed  by  the  cooling  medium,  air  [1 1].  In  equation  (2-33)  an  overall  heat 

watts 

transfer  coefficient,  h  =  30- t  was  used.  This  coefficient  was  derived  for  an  air- 

K  x  meter- 
cooled  machine  and  is  based  on  empirical  calculations  discussed  in  [1 1]. 

Equations  (2-1)  through  (2-32)  are  used  in  Appendix  A  to  do  the  initial  machine 
sizing.  Based  on  the  dimensions  calculated  by  these  equations  and  the  established  machine 
geometry,  the  weight  of  the  machine  was  estimated.  In  addition,  once  the  dimensions  and 
geometry  of  the  machine  have  been  estimated,  other  machine  parameters  such  as  the 
machine  efficiency  and  stator  leakage  reactance  can  be  estimated. 

2.7  Stator  leakage  reactance 

The  stator  leakage  reactance  was  calculated  using  the  methods  of  reference  [11]. 
Leakage  reactance  consists  of  (1)  slot,  (2)  tooth  top,  (3)  end  winding,  and  (4)  harmonic 
reactances.  Tooth  top  leakage  reactance  is  proportional  to  the  gap  length.  In 
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synchronous  machines  with  permanent  magnet  field  the  gap  length  is  small  and  tooth  top 
leakage  reactance  can  be  neglected  [11]. 

For  an  initial  estimate,  the  leakage  reactance  was  assumed  to  consist  of  slot  and 
harmonic  leakage  reactances.  The  slot  leakage  and  harmonic  reactances  were  calculated 
using  the  procedure  outlined  in  [11].  The  calculated  reactance  was  increased  by  10%  to 
account  for  the  effects  of  end  winding  reactance  in  the  total  leakage  reactance. 

The  leakage  reactance  was  assumed  to  be  the  same  for  the  direct  and  quadrature 
axes.  This  assumption  is  commonly  made  in  the  literature  [11]. 

2.8  Losses  and  machine  efficiency 

The  power  dissipated  in  the  machine  determines  its  efficiency.  This  dissipated 
power  will  determine  the  cooling  and  ventilation  requirements  for  the  machine  to  ensure 
that  allowable  temperature  limits  are  not  exceeded.  The  losses  considered  in  Appendix  A 
comprise:  no-load  losses  in  the  iron,  friction  and  windage  losses,  copper  losses  and  load 
losses. 

The  rotating  magnetic  field  will  result  in  heat  losses  in  the  iron  due  to  eddy 
currents  and  hysteresis.  These  type  of  losses  occur  primarily  in  the  stator  teeth,  the 
surface  of  the  rotor  poles  and  the  structural  parts  of  the  machine  exposed  to  alternating 
magnetic  fields.  These  losses  are  proportional  to  the  square  of  the  peak  value  of  B. 
Hysteresis  losses  are  proportional  to  the  frequency  and  eddy  current  losses  are 
proportional  to  the  square  of  the  frequency.  Core  losses  are  estimated  in  Appendix  A 
using  the  relations  discussed  in  [1 1]. 
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Mechanical  losses  in  the  machine  are  the  result  of  friction  in  the  bearings  and 
windage  between  the  rotor  and  the  stator.  These  losses  were  estimated  using  the  relation 
presented  in  [1 1]. 

The  copper  losses  were  calculated  by  estimating  the  resistance  of  the  armature 
windings  and  using  the  rated  terminal  current.  To  calculate  the  armature  resistance,  the 
mean  length  of  conductor  per  phase  in  the  armature  was  estimated  as  suggested  in  [1 1]. 
For  a  three  phase  machine  the  armature  copper  losses  are  3RaI2,  where  Ra  is  the  armature 

winding  resistance  and  I  is  the  line  current. 

Load  losses  are  caused  by  the  flux  produced  by  the  armature  currents.  These 
losses  include  eddy  current  losses  in  the  support  structures,  pole  surfaces  and  damper 
windings.  The  load  losses  were  estimated  as  1%  of  the  armature  copper  losses  [10]. 

The  machine  efficiency  was  estimated  in  Appendix  A  using 

.  losses 

Efficiency  =  1 

input 

where  the  losses  are  the  sum  of  the  individual  losses  described  in  the  preceding 

paragraphs.  The  calculated  machine  efficiency  was  98%,  close  to  the  efficiency  predicted 

by  [6]. 

2.9  Back  iron  sizing 

The  thickness  of  the  back  iron,  he,  must  be  sufficient  to  carry  the  machine  flux. 
Assuming  a  sinusoidal  air  gap  flux  density,  Bgap,  the  back  iron  thickness  is  determined 
using 
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$B-dS  =  0  (2-34) 


If  the  air  gap  flux  is  assumed  to  be  constant,  using  the  specified  machine  geometry 
and  if  R»g,  equation  (2-34)  can  be  simplified  to 


h   = 


gap 


'lO 


(2-35) 


v  Bs  Apy 
where  Bs  represents  the  flux  density  in  the  back  iron. 

The  minimum  back  iron  thickness  is  then  calculated  by  using  a  back  iron  flux 
density  close  to  saturation.  A  value  of  Bs  =  1.2  T,  rms,  was  used  for  the  machine  design. 

2.10  Weight  Calculations 

A  machine  weight  was  estimated  by  calculating  a  rotor  weight  and  a  stator  weight 
for  the  given  machine  geometry.  The  material  composition  of  the  different  components 
was  assumed  to  be  that  of  use  in  standard  motor  construction.  The  machine  weight 
calculated  does  not  include  weight  of  foundations,  the  weight  of  a  cooling  system  or  the 
weight  of  an  enclosure  for  the  motor. 

The  weight  of  the  motor  is  expressed  in  long  tons1  (lton)  for  easier  comparison 
with  current  or  conceptual  drives  for  naval  ships. 


One  long  ton  =  1016  kilograms 
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Chapter  3.  Dynamic  Models 

To  perform  the  fault  simulations  a  model  of  the  permanent  magnet  motor  and  of  the 
electrical  fault  were  derived.  These  two  models  were  incorporated  into  a  dynamic  simulation 
model  developed  by  Professor  James  L.  Kirtley,  MIT.  This  chapter  describes  the  various 
models  used  in  the  simulation.  Models  for  a  permanent  magnet  ac  synchronous  motor,  arc 
fault  and  interconnections  between  the  motor  and  the  power  system  are  presented. 

3.1  Two  axis  transformation 

The  permanent  magnet  motor  model  used  in  this  research  is  based  on  the  derivations 
presented  in  [1 1]  through  [14].  It  assumes  linear  magnetics  and  sinusoidal  stator  winding 
distributions.  Inductance  and  resistance  values  for  the  motor  were  calculated  using  the 
methods  discussed  in  chapter  two  and  Appendix  A. 

Using  the  symmetry  of  cylindrical  rotation,  a  coordinate  transformation  that  accounts 
for  the  relative  motion  between  the  stator  and  the  rotor  is  introduced.  This  transformation 
maps  the  stator  winding  variables  to  a  reference  frame  that  rotates  with  the  rotor.  In  this 
reference  frame,  mutual  inductances  are  independent  of  rotor  position.  Such  transformation  is 
commonly  known  as  the  Park's  transformation.  A  version  of  this  transformation  used  in  [5]  is 
given  by 
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T=    - 


COS0        cos  0 

3 


sin0     -sin  0 


3  +  — 
2tt 
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l 


sin|  0  + 
1 

7? 


(3-1) 


where  0  is  some  arbitrary  angle.  For  a  reference  frame  rotating  with  the  rotor  0  =  cot  where  co 
is  the  speed  of  the  rotor.  This  transformation  is  orthogonal  and  power  invariant, 
regardless  of  the  choice  of  angle.  The  inverse  transformation  is  given  by 


'-i 


T-  =     - 


COS0 


cos  0 

cos  0  +  — 


sin0 


sin  0-- 


sin  0  + 


2tc 

3 

2tt 
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(3-2) 


If  the  Park's  transformation  is  applied  to  an  arbitrary  vector,  F,  representing  a  set  of 
stator  variables,  such  as  currents  or  voltages,  the  stator  variables  are  mapped  to  a  fixed  rotor 
reference  frame.  In  this  stationary  reference  frame  the  machine  reactances  are  constant, 
independent  of  rotor  position. 


^dqO  ^abc 

F     =  T  'F 

1  abc  1       A  dqO 


(3-3) 


The  d  and  q  subscripts  are  used  to  represent  the  direct  and  quadrature  axes, 
respectively.  The  direct  axis  is  aligned  with  the  polar  axis  and  the  quadrature  axis  is  aligned 
with  the  interpolar  space.  Zero-sequence  variables  are  represented  by  the  subscript  0  in 
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equations  (3-3).  The  zero-sequence  components  represent  components  of  armature  currents 
that  do  not  produce  a  net  air  gap  flux  [4], 

3.2  Per  Unit  Scaling 

The  models  used  in  the  simulation  have  been  scaled  to  a  per  unit  (p.u.)  system.  Scaling 
to  a  per  unit  system  produces  variables  whose  magnitudes  are  close  to  one.  The  base  values 
selected  for  this  scaling  can  be  arbitrary.  However,  for  this  analysis  it  is  convenient  to  use  the 
motor's  rated  power  and  voltage  as  the  base  values  for  the  per  unit  scaling.  Expressing  the 
motor  variables  in  a  per  unit  system  allows  the  comparison  of  this  motor  against  other  motors 
regardless  of  their  ratings. 

For  this  motor  the  base  quantities  are  defined  by 

P   =— V  I 

rB         «    VB1B 

ZB=^  (3-4) 

*B 

A,B  =  — 

To  convert  the  motor  parameters  to  the  per  unit  system,  the  ordinary  variables  were 
divided  by  the  corresponding  base  quantities.  A  new  quantity  will  be  introduced  into  the  model 
to  change  torque  to  the  per  unit  system.  This  new  constant  is  defined  as  the  inertia  constant, 
H;  it  represents  the  rotor's  kinetic  energy  at  rated  speed  divided  by  the  motor  rated  power. 
The  inertia  constant,  H,  has  units  of  seconds  and  is  expressed  as 

1 


Jco 
P 


H  =  ^—  (3-5) 


B 
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In  equation  (3-5),  J  is  the  moment  of  inertia  of  the  shaft  and  com  is  the  shaft's  rotational 
speed.  In  general  the  moment  of  inertia  of  a  circular  shaft  can  be  calculated  as  J  =  f  rpdV , 

V 

where  p  is  the  shaft's  material  density  [3].  It  can  be  shown  that  for  a  circular  shaft  this  is 
equivalent  to  J  =  Wk: ,  where  W  is  the  mass  of  the  shaft  and  k  is  its  radius  of  gyration.  For  a 
solid  circular  shaft  of  radius,  R,  the  radius  of  gyration  is  given  by  k  =  R./V2  .  A  value  for  H 
was  calculated  for  this  motor  using  the  mass  and  rotor  radius  calculated  in  Appendix  A. 

Since  the  rotor  is  directly  connected  to  the  propulsion  shaft,  the  length  of  the  shaft  and 
the  propeller  needs  to  be  included  in  the  moment  of  inertia  calculation.  The  weight  of  the 
propeller  and  of  the  propulsion  shafts  can  be  estimated  from  propeller  and  shaft  weights  from 
U.S.  Navy  ships  or  commercial  ships. 

3.3  Permanent  Magnet  Motor  Model 

The  permanent  magnet  machine  model  is  defined  by  a  set  of  six  coupled  electro- 
mechanical equations.  The  electrical  equations  in  this  set  are  statements  of  Kirchoff  and 
Faraday's  laws,  which  describe  the  voltage-current  relations  in  the  machine.  The  mechanical 
equations  are  statements  of  Newton's  Second  Law  of  Motion.  Coupling  between  the 
mechanical  and  electrical  systems  is  through  the  dependence  of  electrical  torque  on  flux 
linkages  and  through  the  dependence  of  flux  linkages  on  torque  angle,  5. 

The  machine  parameters  used  for  the  simulations  are  calculated  in  Appendix  A  using 
the  methodology  discussed  in  chapter  two.  A  two-axis  (d-q  axis)  representation  of  the 
synchronous  machine  is  used. 
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In  the  d-q  reference  frame,  the  direct  axis  of  the  motor  is  represented  by  the  following 
circuit  based  on  the  model  presented  in  references  [5]  and  [12]. 

Ra  Xal 

-A\ vQGQ/ — 


+ 


Vd 


Id 


Xkd 


Xad 


(V\   If  n 


Rkd 


Figure  3-1.  Direct  axis  circuit  representation 

where  the  constant  current  source  represents  the  field  generated  by  the  permanent  magnets. 
Similarly  the  quadrature  axis  of  the  motor  can  be  represented  by  the  following  circuit. 

Ra  Xal 
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Figure  3-2.  Quadrature  axis  circuit  representation 

The  following  notation  has  been  introduced  in  the  above  two  models: 

Xad  and  Xaq  represent  the  magnetizing  reactances, 

Xai  represents  the  armature  leakage  reactance,  which  is  assumed  to  be  equal 
for  the  d  and  q  axes, 

Ra  is  the  armature  resistance, 
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Xkd  and  Xkq  are  the  damper  winding  reactances,  and 

Rkd  and  Rkq  are  the  damper  winding  resistances. 

Using  the  models  shown  in  figures  3-2  and  3-3,  the  three  phase  ac  synchronous  machine  is 
represented  by  the  following  set  of  equations 

V,  = — ^  +  R  l.-X  co 

d  ^  a     d  q 

dX 

vq=-^-+Raiq+V> 

Vkd=^  +  RkdIkd  (3-6) 

V     —  q  4-R     T 

Vkq    -        j  "1"IVkq1kq 

t.=^(v,-V*)— 

2  a>_ 

m 

For  permanent  magnet  excitation,  the  field  is  represented  by  the  constant  current  source ,  Itm- 
The  flux  equations  for  the  machine  are 


^d=LdId+LadIkd+LadIrm 
^kd   =  Lad^d  +^kd^kd  +^ad^fm 


\=LqIq+LaJkq 
^kq    =LaqIq   +  LkqIkq 


where  the  inductances,  Lj  and  Lq,  are  defined  by 


(3-7) 


Ld  =L„  +Lad 


Lq  =Lal+Laq 


The  simulation  code  was  written  to  be  used  with  any  ac  synchronous  machine, 
regardless  of  its  rating.  For  this  reason  it  is  intended  to  be  used  with  machines  scaled  in  the  per 
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unit  system.  To  use  the  dynamic  simulation  code,  the  above  set  of  equations  were  scaled  in  the 
per  unit  system.  Dividing  by  the  appropriate  base  quantities,  the  sets  of  equations  (3-5)  and  (3- 
6)  are  represented  in  the  per  unit  system  by 


Vd 

Vkd 

M>kq 


id 

xd 

Xad 

Xad 

Xad 

Xkd 

Xad_ 

!kd 

X,  Xaq 

Xaq         Xkq 


Lkq 


(3-8) 


(3-9) 


1    dvj/ 


CD 

+  ra1d-Vq  — 
co  „    dt  co  „ 


1    dvj/c 


A,       +I"a1q+V|/d 

co„    dt 


co 
co„ 


1   d\|/ 


kd 


co„    dt 


0 


1    dv|/ 


kq 


+  I'kd1kd 


+  r,,  i, 


(3-10) 


co0     dt        kqkq 
Te=(M>diq-Myd) 

The  damper  windings'  voltages,  v^d  and  v^q,  were  set  to  zero  since  these  windings  are 
normally  shorted.  Using  a  Thevenin  equivalent  circuit  of  the  d-axis  model,  figure  3-1,  the 
constant  current  source  can  be  represented  as  a  voltage  source,  Eal  =  XadIfm  in  series  with  the 

magnetizing  impedance,  Xa(j. 
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Figure  3-3.  Thevenin  equivalent  of  d-axis  model 

With  the  model  presented  in  figure  3-3,  equation  (3-7)  can  be  rewritten  as  follows 

(3-11) 

By  solving  the  above  equations  for  the  currents,  the  d-q  model  of  the  motor  can  then 
be  incorporated  into  a  simulation  code  that  connects  the  motor  to  an  external  network. 
Solving  equations  (3-9)  through  (3-11)  for  the  currents,  id  and  iq,  and  after  some  simplifications 
the  following  sixth  order  state-space  representation  of  the  motor  was  obtained 


dt 


V< 


+  -5-  +  G>\|/q+G>0Vd        (3-12) 


T        T 

xad  xad 


dM/q         v|/( 


—  =  -— ^-coi|/d+co  v        (3-13) 

dt  T        T  q 

aq  aq 


de"  e"      (x    -x")         e 

q   _  q        \Kd       Ad/;      ,    caf 

1  A        ' 


dt 


de^ 
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t  n  -r  n 

1do  Xdo 
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rrri 
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d5 

dt 


=  CO  -co 


(3-14) 


(3-15) 


(3-16) 
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=  — -(Tm-T)  (3-17) 

dt      2H     m       e7  ' 


In  equations  (3-12)  through  (3-17)  the  following  quantities  are  defined: 


e"  =  — —  vj/kd  =  Voltage  behind  subtransient  reactance, 

xkd 


e"  ~  — ~  H'kq  =  Voltage  behind  subtransient  reactance, 
x 

q 


Id"  =  — ~  =  D-axis  open  circuit  subtransient  time  constant, 

°>orkd 


To"  =  — ~ =  Q-axis  open  circuit  subtransient  time  constant, 

Q  r,  . 

o    kd 


X 

Tad  =  — —  =  Direct  axis  armature  time  constant, 


K 
T    =  — —  =  Quadrature  axis  time  constant, 

Qora 


xd'  =  xd —  =  D-axis  subtransient  reactance, 


kd 


Xaq 

x"  =  xn =  Q-axis  subtransient  reactance, 

x 

Xkq 


The  stator  currents,  in  the  motor  reference  frame,  have  been  defined  by 

V.  -e" 

id=^V  <3"18) 

xd 

\\j   +  e" 
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Equations  (3-12)  and  (3-13)  are  the  stator  equations,  which  have  time  constants  of  the 
order  of  1/co  0  =  0.0026  seconds.  If  the  transients  of  interest  are  in  the  order  of  0. 1  to  several 
seconds  then  two  stator  equations,  equations  (3-12)  and  (3-13),  will  become  algebraic 
equations  provided  that  the  following  conditions  hold 

1       1      d 


co  » 


Tad      Taq      dt 


Furthermore  if  co  *  co  o  then  equations  (3-12)  and  (3-13)  simplify  to 

Vd=vd  (3-12a) 

Vq=vq  (3-13a) 

With  these  two  algebraic  equations  the  machine  model  has  been  reduced  to  a  fourth  order 
model,  equations  (3-14)  through  (3-17). 

3.4  Network  model 

The  external  network  model  used  in  the  simulation  code  is  based  on  the  model  derived 
in  [15].  This  model  uses  the  machine  currents  to  interface  with  the  network.  The  definition  of 
the  network  is  accomplished  by  defining  the  nodes  and  branches  of  the  network.  This  model  is 
discussed  in  chapter  four. 

3.5  Fault  Model 

The  fault  model  used  for  this  research  is  an  ac  arc  fault  whose  properties  and 
characteristics  are  described  in  detail  in  [2].  Reference  [2]  describes  the  electric  arc  as  a  self- 
sustained  electrical  discharge  having  a  low  voltage  drop  and  capable  of  supporting  large 
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currents.  At  atmospheric  pressure  the  temperature  of  the  arc  will  reach  temperatures  as  high  as 
6,000K. 

According  to  [2],  for  large  currents  the  arc  voltage,  ea,  is  given  by 


eQ  =  c,i  +-7-  +  C3I' 


C2 

i 


(3-20) 


where  the  constants  c\  through  C3  are  determined  by  the  electrode  materials.  The  voltage- 
ampere  characteristic  of  ac  arcs  will  show  hysteresis  effects  with  distinct  ignition  and  extinction 
voltages  [2]. 

The  arc  is  modeled  as  a  constant  voltage  drop;  for  calculation  simplifications  the  values 
of  ignition  and  extinction  voltages  will  be  neglected.  For  the  materials  considered  such  as 
copper  and  iron  the  arc  voltage  drop  is  approximately  60  volts  [2].  The  voltage-current 
relation  of  the  fault  can  then  be  approximated  by  a  characteristic  such  as  that  shown  in 

figure  3-4. 

v 


ea 


Figure  3-4.  Fault  voltage-current  characteristic 


40 


For  the  purposes  of  this  research  the  impedance  of  the  fault  is  assumed  to  be  much 
smaller  than  the  synchronous  impedance  of  the  motor.  This  is  simulated  by  using  fault 
impedance  values  ten  times  lower  than  the  synchronous  impedance  values  for  the  motor. 

The  power  dissipated  by  the  fault  is  estimated  by  multiplying  the  fault  current  by  the 
fault  voltage. 

3.6  Motor  Load  Model 

The  mechanical  load  placed  on  the  motor  by  the  ship  is  represented  by  a  mechanical 
torque  on  the  shaft.  The  power  delivered  by  a  propulsion  shaft  for  a  ship  moving  at  high  speed 
is  proportional  to  the  speed  of  the  ship  cubed,  P  x  v3  [16].  This  assumption  is  not  valid  for 
ships  such  as  destroyers  that  can  reach  speeds  near  or  above  hull  speed2.  At  speeds  close  to 
the  hull  speed,  the  power  required  may  vary  with  speed  raised  to  a  power  approaching  6  or  7 
[16]. 

Once  a  propeller  has  been  selected  and  matched  to  the  hull,  it  can  be  shown  the 
rotational  speed  of  the  shaft  is  proportional  to  the  speed  of  the  ship.  Using  this  relation,  and 
since  P  =  Tco  m  then  the  torque  on  the  shaft  can  be  assumed  to  be  proportional  to  the  square 

of  the  shaft  speed.  This  torque  is  assumed  to  be  of  the  form 

T   =aor  (3-21) 

mm  v  • 


m  m 


2  Hull  speed  or  critical  speed  length  ratio  is  the  speed  at  which  wave  making  resistance  starts 
accumulating  most  rapidly.  Hull  speed  is  calculated  as  Vhull  (knots)  =  1.34^/E^  ,  where  Ls  is  the  ship's 

length  in  feet. 
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The  constant  is  calculated  with  the  motor  operating  at  full  load  using  the  rated  shaft 
speed  and  torque  calculated  in  Appendix  A.  It  is  further  assumed  that  this  constant  is 
independent  of  shaft  speed. 

A  windmilling  shaft  is  not  delivering  any  thrust  or  torque;  however,  it  will  rotate  at  a 
speed  determined  by  the  characteristics  of  the  propeller  and  the  speed  of  the  ship.  Once  the 
propeller  is  selected,  for  a  fixed  pitch  propeller,  the  speed  of  rotation  of  the  shaft,  while 
windmilling,  will  be  proportional  to  the  speed  of  the  ship.  So  as  the  ship  slows  down  the 
windmilling  shaft  will  slow  down  proportionately.  To  incorporate  this  model  into  the  dynamic 
simulation  model  it  is  necessary  to  determine  how  the  ship's  speed  changes  as  a  function  of 
time  while  the  ship  is  coasting  down. 

For  a  ship  moving  on  a  constant  course  at  some  initial  speed,  v0,  the  total  kinetic 

energy  of  the  ship,  U,  is  of  the  form  U  =  —  Mv;: .  The  symbol  M  accounts  for  the  mass  of  the 

ship  and  the  "added"  mass  in  the  direction  of  motion.  While  a  ship  coasts  down  from  some 
given  speed  the  rate  at  which  the  ship  loses  energy  to  the  water  is  equal  to  the  effective 
horsepower,  EHP,  of  the  ship 

—  =  -EHP  (3-22) 

dt 

In  general  EHP  is  proportional  to  speed  elevated  to  some  power,  EHP  =  kvn ,  where  k  is 

some  constant.  During  a  short  time  interval,  At  =  t„  - 1 ,  equation  (3-22)  can  be  written  after 

some  manipulations  as 

-M(v;-v2)  =  -(t()-t)(kvn)  (3-23) 


42 


If  the  ship  is  assumed  to  be  moving  at  an  initial  speed,  v0,  at  time,  t0  =  0,  then  from  (3-23) 

a       2kvn 
v;-v=— -t  (3-24) 

M 

For  speeds  close  to  the  initial  ship  speed,  v  »  v0 ,  equation  (3-24)  can  be  rewritten  as 

*-*£!,  ,3-25, 

M 

Using  equation  (3-25)  the  following  relation  between  ship's  speed  and  time  is  obtained 

v  x  — -  (3-26) 

n^t  V        ' 

Since  the  rotational  speed  of  a  windmilling  shaft  is  proportional  to  the  speed  of  the 
ship,  using  equation  (3-26)  the  rotational  speed  of  the  shaft,  com,  will  exhibit  the  same  time 
dependence  as  the  ship's  speed. 

If  the  ship  is  operating  below  hull  speed,  which  is  the  case  for  most  full  form  ships,  then 

at  its  rated  speed  it  can  be  assumed  that  EHP  x  v3 .  For  n=3,  then  v  oc  — ■=.  After  the 

Vt 

casualty  takes  place  and  the  motor  main  supply  breakers  are  opened,  the  shaft  will  continue  to 

rotate  at  a  speed  proportional  to  the  speed  of  the  ship.  It  can  be  assumed  that  the  increase  in 

ship's  resistance  added  by  the  windmilling  shaft  is  negligible  compared  to  the  hull  resistance 

[17]. 

Considering  the  discussion  above,  while  the  ship  is  coasting  down,  the  rotational  speed 

of  the  shaft  was  assumed  to  be  of  the  form  co  m  x  —= .  The  final  shaft  speed  will  be 

v  t 

proportional  to  the  final  ship  speed. 
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Chapter  4.  Simulation  Model 

With  the  system  model  complete,  as  discussed  in  chapter  three,  the  dynamic 
response  of  the  permanent  magnet  motor  and  the  internal  fault  was  investigated.  The 
purpose  of  conducting  these  studies  was  to  determine  the  effect  of  the  internal  fault  after 
the  motor  is  disconnected  from  the  system.  The  possible  damage  done  by  the  fault  will 
occur  in  a  very  short  time  so  that  the  time  delay  of  the  relay  and  protection  system  might 
not  be  significant. 

A  possible  way  to  reduce  the  effects  of  the  fault  is  to  short  circuit  the  motor 
terminals  after  the  motor  is  disconnected  from  its  power  supply.  By  shorting  the  motor 
terminals  it  might  be  possible  to  reduce  the  power  input  into  the  fault,  thus  serving  to 
minimize  the  effects  of  the  fault. 

4.1  Discussion  of  the  simulation  program 

The  simulation  code  was  written  in  C++.  Using  this  language  allowed  the  use  of 
object-oriented  programming  (OOP).  This  type  of  programming  consists  of  building 
programs  as  a  collection  of  abstract  data  type  instances.  Operations  performed  on  the 
object  types  are  the  abstract  operations  that  solve  the  problem.  These  objects  serve  as 
modules  that  can  be  reused  for  solving  another  problem  in  the  same  domain  [17].  In  the 
simulation  code  generators,  motors,  nodes  and  lines  are  built  into  objects.  These  objects 
are  invoked  throughout  the  code  as  needed  to  solve  the  simulation  problem. 
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The  dynamic  simulation  program,  Appendices  B  through  E,  requires  the  user  to 
provide  four  input  files.  The  first  two  files  are  used  to  define  the  node  types  and  the 
interconnection  lines  between  nodes  and  their  reactances.  In  addition,  they  define  the 
nodes  to  which  the  generators  or  motors  are  connected  and  the  names  of  the  output  files 
that  will  contain  the  node  voltages  and  line  currents. 

Two  other  input  files  are  needed  for  the  simulation.  One  file  defines  the  electrical 
parameters  of  the  machine.  The  other  file  defines  the  timing  and  sequence  of  events  for 
the  simulation.  The  machine's  electrical  parameters  used  in  the  simulation  were  those 
calculated  using  Appendix  A. 

The  simulation  code  consists  of  four  executable  programs  that  can  be  run 
individually  or  using  a  shell  program  that  executes  all  four  programs  in  their  appropriate 
sequence.  The  first  program  executed  is  a  pre-processor  program.  It  formats  the  node 
and  line  input  data  into  a  format  that  can  be  used  by  the  rest  of  the  code.  The  next 
program  in  the  simulation  estimates  the  power  flow  in  the  network.  It  assigns  initial 
values  to  node  voltages  and  line  currents.  Using  the  power  flow  equations,  the  voltages 
and  current  flows  of  the  connected  system  are  calculated. 

The  third  executable  program  uses  the  node  voltage,  line  current  and  power  flow 
data  to  build  the  line  data  required  by  the  simulation  program.  The  simulation  program 
will  generate  three  output  files  that  contain  the  state  variables,  the  node  voltages  and  the 
line  currents.  To  perform  the  simulation  the  four  programs  need  to  be  run  in  sequence 
since  the  output  files  of  one  program  serve  as  input  files  for  the  next  program. 
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The  dynamic  simulation  program,  developed  by  Professor  James  L.  Kirtley,  MIT, 
was  originally  written  to  simulate  the  dynamic  response  of  wound  field  ac  synchronous 
generators.  The  generators  are  defined  as  objects  in  the  code  so  that  multiple  generators, 
all  possibly  different  can  be  connected  to  the  network.  For  this  research  it  was  necessary 
to  generate  a  similar  object  for  the  permanent  magnet  motor.  To  properly  interface  the 
motor  object  with  the  other  programs  it  was  necessary  to  maintain  the  same  notation  as 
that  used  for  the  generator  objects. 

The  simulation  programs  were  modified  to  incorporate  the  permanent  magnet 
motor.  To  incorporate  the  motor  into  the  simulation  program  the  following  major 
modifications  were  made  to  the  simulation  code: 

1.  The  state  equations  for  the  machine  were  modified  to  incorporate  the 
permanent  magnet  motor  model  discussed  in  chapter  three.  The  reference  frame  of  the 
model  was  changed  to  the  motor  reference  frame. 

2.  The  programs  were  modified  to  delete  the  voltage  regulator  associated  with 
each  generator. 

3.  The  ship  model  discussed  in  chapter  three  was  incorporated  into  the  program. 
This  model  is  used  to  drive  the  windmilling  shaft  after  the  motor  is  disconnected  from  its 
power  supply. 

4.  The  program  was  modified  to  allow  connecting  a  motor  to  the  network. 

5.  The  fault  model  was  incorporated  into  the  motor  simulation  code.  The  fault  is 
introduced  at  a  predetermined  time  after  the  simulation  starts. 

6.  The  simulation  code  was  modified  to  output  fault  current  as  a  function  of  time. 
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4.2  Fault  simulation 

The  permanent  magnet  motor  model  and  the  parameters  derived  in  chapters  two 
and  three  were  incorporated  into  the  dynamic  simulation  model,  Appendix  E.  The 
simplified  system  shown  in  figure  4-1  was  used  for  the  simulations.  For  this  simulation  it 
was  assumed  that  a  single  permanent  magnet  motor  was  connected  to  an  infinite  bus.  The 
idea  of  an  infinite  bus  is  not  applicable  to  ships'  power  plants  because  loads  such  as 
propulsion  motors  can  have  ratings  comparable  to  that  of  the  generators. 

The  infinite  bus  assumption  was  made  to  establish  the  initial  state  of  the  system. 
For  the  simulation,  the  motor  is  operating  at  rated  power  and  rated  speed.  At  the  time  of 
the  fault  the  affected  motor  was  disconnected  from  the  rest  of  the  system  so  that  there  was 
no  interaction  between  the  motor  and  the  rest  of  the  system.  Since  interaction  effects  will 
not  be  addressed,  the  concept  of  an  infinite  bus  was  considered  suitable  for  establishing 
the  initial  state  of  the  system.  This  bus  represented  the  ship's  generators  whose  output  can 
be  adjusted  to  provide  a  specific  voltage  and  power  flow. 

Using  the  model  shown  in  figure  4-1,  the  motor  was  disconnected  from  the  system 
some  time  after  the  initiation  of  the  fault.  In  a  real  system  this  would  have  been 
accomplished  with  circuit  breakers  together  with  some  sort  of  sensing  system.  For  this 
simulation  a  fixed  time  of  100  milliseconds,  after  the  initiation  of  the  fault,  was  chosen  for 
the  breakers  to  trip  and  disconnect  the  motor  from  its  power  supply.  This  time  was  used 
for  all  simulation  studies. 
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Once  the  motor  was  disconnected  from  its  power  supply,  it  was  allowed  to 
windmill  as  the  shaft  coasted  down.  The  speed  of  the  windmilling  shaft  was  taken  to  be  a 
function  of  the  speed  of  the  ship,  as  previously  discussed  in  chapter  three. 

The  fault  current  and  power  dissipated  at  the  fault  were  calculated.  The 
simulations  were  terminated  after  one  second.  This  time  was  selected  for  the  simulations 
since  as  will  be  shown  most  of  the  power  dissipated  at  the  fault  will  occur  within  this  time. 
This  is  consistent  with  the  results  obtained  in  reference  [19].  In  addition,  this  short  time  is 
consistent  with  the  assumptions  made  in  chapter  three  for  the  ship  model. 
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Figure  4-1.  System  model 

The  internal  fault  in  the  motor  was  simulated  by  assuming  that  the  fault  occurred 
somewhere  in  the  windings  between  two  points  which  have  generated  internal  voltages 
and  leakage  inductances.  A  similar  model  was  used  to  simulate  internal  faults  in  [19]. 

The  motor  and  power  supply  interconnection  model  is  shown  in  figure  4-2.  This 
model  consists  of  an  internal  voltage  source  and  in  series  with  a  reactance. 
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Figure  4-2.  Motor  and  Network  Interconnection  Model 

A  derivation  for  this  model  and  its  representative  equations  are  discussed  in  [15]. 
For  this  model  it  is  assumed  that  subtransient  saliency  can  be  neglected,  x^  =  x" .  The 

internal  voltages,  ea",eb",ec"  ,  represent  voltages  induced  by  rotor  fluxes  and  xg 

represents  the  impedance  in  the  neutral  of  the  machine  and  of  mutual  coupling  between  the 
phases.  The  internal  voltages  are  calculated  as  follows 


e   =- 


e    = 


e  = 


—  (e"  sinG-  e?cos0)  +  — cose^s.  +  —  sinG  ^L 
o0  co0  dt      (o„  dt 

(e/'sin(0-— )-e^'cos(0-— ))  + — cos(0-— )— ^  +  - -sin(0-— )— ^ 

coo  3  3         coo  3      dt      coo  3dt 

co       ,  2tt.       „  2tt.       1         ..     27C.de;'       l  2tt  de? 

—  (e"  sin(0  +  — )  -  o,"A  cos(0  +  — ))  +  —  cos(0  +  — )  — -  +  —  sin(0  +  — )  — - 
co  „  3  3         co ,  3      dt      co  „  3      dt 


where,  0  is  the  shaft  angle  and  e  d  and  e  q  are  the  subtransient  voltages  defined  in  chapter 
three.  A  detailed  derivation  of  this  model  is  explained  in  [12]  and  [15]. 

This  machine  model  was  used  in  the  synchronous  machine  simulation  program. 
The  internal  fault  was  assumed  to  occur  between  two  phase  points  whose  voltages  are 
proportional  to  the  internal  phase  voltage.  Because  the  fault  occurred  between  two 
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phases,  the  voltage  across  the  fault  is  proportional  to  the  difference  between  the  internal 
voltages  of  the  two  points. 

The  internal  voltages  in  the  machine  are  generated  across  the  entire  winding. 
Since  the  fault  can  occur  anywhere  in  the  windings  there  are  many  possible  connections 
that  can  be  established.  For  this  analysis  the  fault  was  assumed  to  occur  in  the  middle  of 
the  winding.  The  impedance  seen  by  the  fault  will  be  proportional  to  the  impedance  of  the 
motor  winding. 

The  power  dissipated  in  the  fault  was  calculated  assuming  that  the  voltage  drop 
across  the  fault  was  constant  as  discussed  in  chapter  three. 

This  type  of  fault  and  simulations  for  a  large  turbogenerator  are  discussed  in  [19]. 
Fault  currents  as  large  as  14  per  unit  and  peak  phase  currents  of  approximately  7  per  unit 
were  obtained  in  [19]  for  an  internal  phase  to  phase  fault.  These  large  currents  can  be 
expected  to  cause  severe  stator  damage. 

For  an  asymmetrical  fault,  large  negative  phase  sequence  currents  are  expected  to 
flow  in  the  stator  windings  and  consequently  in  the  rotor  damper  bars.  These  large 
currents  in  the  damper  bars  could  impose  high  thermal  stresses  in  the  damper  bars. 

4.3  Simulations 

The  following  cases  were  run  for  this  research  and  the  results  obtained  are 
presented  and  discussed  in  chapter  five.  All  the  simulations  were  started  with  the  motor 
operating  at  rated  speed  and  power.  A  phase  to  phase  fault  was  introduced,  in  all  cases, 
300  milliseconds  after  the  simulation  was  initiated.  The  motor  was  disconnected  from  its 
power  supply  0.1  seconds  after  the  initiation  of  the  fault. 
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For  the  first  simulation  the  motor  terminals  remained  open  after  the  motor  was 
disconnected  from  its  power  supply.  Once  the  motor  was  disconnected  it  was  allowed  to 
windmill  as  the  ship  coasted  down,  as  described  in  chapter  three. 

In  the  next  set  of  simulations  the  motor  terminals  were  shorted  as  soon  as  the 
motor  was  disconnected  from  its  power  supply.  The  motor  was  allowed  to  coast  down 
with  the  terminals  shorted  for  the  remaining  time  of  the  simulation. 

The  same  two  cases  discussed  above  were  run  for  25%  and  75%  of  the  winding 
extents. 

4.4  Fault  Parameters 

Consistent  with  the  fault  model  described  in  chapter  three  the  voltage  across  the 
fault,  ea,  was  maintained  constant  at  a  value  of  0.05  per  unit.  This  value  is  based  on  the 

arc  voltage  values  for  copper  electrodes  given  in  [2]  and  the  base  voltage  of  Appendix  A. 
The  fault  currents  calculated  in  the  model  are  consistent  with  internal  fault  currents  as 
discussed  in  [19]. 

The  fault  current  is  calculated  using  the  model  presented  in  [2].  For  a  constant 
fault  voltage  the  relation  between  fault  current  and  voltage  is  given  by 

L—  +  Ri  +  e    =v  (4-1) 

dt 

where  v  is  the  applied  voltage.  Equation  (4-1)  represents  a  series  R-L  circuit  with  an 
applied  voltage,  v,  and  a  constant  voltage  drop,  ea,  that  represents  the  burning  voltage  of 
the  arc.  For  these  simulations  the  applied  voltage  is  proportional  to  the  difference 
between  two  winding  voltages. 
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If  the  resistance,  R,  is  neglected  compared  to  the  inductive  reactance  of  the  circuit, 
and  assuming  ea  is  constant,  then  equation  (4-1)  is  solved  for  the  fault  current,  isc,  [2] 

isc  = -cos(cot +  <}>)  +  —  (--tot) 

toL  toL   2 

where,  (b  =  cos"1  — — . 

m 

The  instantaneous  power,  psc,  dissipated  by  the  fault  was  calculated  using  the 
product  of  fault  current  times  fault  voltage,  or  psc  =  eaisc . 
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Chapter  5.  Conclusions  and  Results 

5.1  Simulation  Results 

The  following  sections  describe  the  simulations  that  were  run  for  this  research. 
For  the  first  set  of  simulations  the  motor  was  disconnected  from  its  power  supply  1 00 
milliseconds  after  the  fault  was  initiated.  The  terminals  of  the  motor  remained  open 
circuited  as  the  motor  windmilled.  Three  cases  were  simulated  with  the  fault  occurring  in 
different  locations  of  the  winding. 

The  same  cases  discussed  above  were  run  with  the  terminals  of  the  motor  shorted 
at  the  time  the  motor  supply  breakers  opened.  It  was  assumed  that  the  motor  terminals 
were  shorted  at  the  same  instant  that  the  motor  was  disconnected  from  its  power  supply. 

5.1.1  Simulations  with  motor  terminals  open 

For  these  simulations  the  motor  terminals  remained  open  after  the  supply  breakers 
to  the  motor  were  opened.  The  internal  fault  occurred  300  milliseconds  after  the 
simulation  started.  The  motor  was  disconnected  from  its  power  supply  1 00  milliseconds 
later  and  allowed  to  windmill.  Results  from  this  simulation  are  shown  in  figures  5-1 
through  5-8.  These  figures  correspond  to  a  fault  occurring  in  the  middle  of  the  windings. 

The  arc  current  reached  peak  levels  of  approximately  5  per  unit,  figure  5-1,  while 
the  motor  was  connected  to  the  power  supply.  After  the  motor  was  disconnected  the  fault 
current  decreased  to  peak  values  of  nearly  2  per  unit.  As  the  motor  slows  down  the 
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internal  voltage  generated  inside  the  motor  will  decrease,  equation  ( 1  - 1 ).  A  decreasing 
internal  voltage  results  in  decreasing  fault  currents,  equation  (4-1). 

The  arc  voltage  is  shown  in  figure  5-2.  This  voltage  corresponds  to  a  constant  arc 
burning  voltage.  Using  the  method  described  in  chapter  four  the  instantaneous  power 
dissipated  at  the  fault  was  calculated,  figure  5-3.  For  this  fault  the  average  power 
dissipated  during  the  transient  was  determined  to  be  approximately  0.06  per  unit. 
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Figure  5-1.  Arc  current  as  a  function  of  time 
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Figure  5-2.  Arc  voltage  as  a  function  of  time 
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Figure  5-3.  Power  dissipated  in  arc 


Figure  5-4.  Voltage  behind  subtransient  reactance 


Figure  5-5.  Voltage  behind  subtransient  reactance 

Figures  5-6  through  5-8  show  in  more  detail  the  characteristics  of  the  arc  voltage, 

current  and  instantaneous  power  dissipated  at  the  fault. 
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Figure  5-6.  Arc  voltage 
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Figure  5-7.  Arc  current 
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Figure  5-8.  Arc  power 

The  arc  current,  voltage  and  power,  shown  in  figures  5-9  through  5-11,  are  for  a 

location  in  the  machine  where  the  peak  voltage  is  one  fourth,  25%  winding  location,  of  the 
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internal  voltage.  The  average  power  dissipated  in  this  case  over  the  duration  of  the 
transient,  700  milliseconds,  was  0.035  per  unit  with  a  maximum  current  peak  of  3  per  unit. 


73 

ci, 

c      0 

>-< 

5 


f\l\l\l\ 


a   a  A  A  A  A  A  A  A  A  A  A 

/  t/  1/ 


v  V  li   1/  U    v 


0.3 


0.35 


0.4 
Time, sec 


V    V 


0.5 


Figure  5-9.  Arc  current 
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Figure  5-10.  Arc  power 
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Figure  5-11.  Arc  voltage 
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Figures  5-12  through  5-14  correspond  to  a  fault  at  a  location  where  the  maximum 
voltage  is  three  fourths,  75%  winding  location,  of  the  generated  internal  voltage.  For  this 
case  the  average  power  dissipated  at  the  arc  was  0.09  per  unit  with  a  maximum  peak 
current  of  6.2  per  unit. 


Figure  5-12.  Arc  voltage 
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Figure  5-13.  Arc  power 
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Figure  5-14.  Arc  current 
5.1.2  Simulations  with  shorting  of  the  motor  terminals 

The  purpose  of  these  simulations  was  to  determine  if  shorting  the  motor  terminals, 
after  the  motor  has  been  electrically  disconnected,  has  any  effect  on  the  power  dissipation 
in  the  fault.  Figures  5-15  through  5-22  show  the  arc  currents,  voltages  and  power 
dissipated  at  the  fault.  The  terminals  of  the  motor  were  shorted  immediately  after  the 
motor  supply  breakers  were  opened. 


Figures  5-15  through  5-17  correspond  to  fault  locations  near  the  middle  of  the 


winding. 
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Figure  5-15.  Arc  current  as  a  function  of  time 
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Figure  5-16.  Arc  voltage  as  a  function  of  time 


Figure  5-17.  Power  dissipated  at  fault 

Shorting  the  terminals  of  the  motor  reduced  the  average  power  dissipated  in  the 
arc  to  0.051  per  unit  with  a  peak  current  of  4.9  per  unit.  This  represents  an  approximate 
20%  reduction  in  the  power  dissipated  in  the  arc. 

Figure  5-18  through  5-20  correspond  to  a  25%  winding  location.  For  a  fault  at 
this  location  the  average  power  dissipated  was  0.0295  per  unit  with  a  peak  current  of  3 
per  unit. 
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Figure  5-18.  Arc  current 
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Figure  5-19.  Arc  voltage 
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Figure  5-20.  Arc  power 
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For  a  75%  winding  location  the  simulation  results  are  shown  in  figures  5-21 
through  5-23.  In  this  case  the  average  power  dissipated  was  0.08  per  unit  with  a  peak 
current  of  6.2  per  unit 
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Figure  5-21.  Arc  current 
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Figure  5-22.  Arc  power 
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Figure  5-23.  Arc  Voltage 

The  simulations  where  the  motor  terminals  were  shorted  showed  a  decrease  in  the 

average  power  dissipated  at  the  fault.  For  all  three  cases,  an  approximate  20%  reduction 
in  dissipated  power  was  observed. 


5.2  Limitations  of  the  simulation  models 

The  models  used  in  this  research  did  not  address  interactions  between  other 
generators  and  loads  during  the  time  of  the  fault.  In  naval  propulsion  plants,  the  size  of 
the  generators  can  be  of  the  same  magnitude  as  some  of  the  electrical  loads.  The  effects 
of  an  internal  fault  on  the  rest  of  the  electrical  system  should  be  considered.  Such  studies 
could  be  accomplished  by  incorporating  the  models  developed  in  this  research  into 
simulation  programs  such  as  those  in  [4]. 

The  only  case  considered  in  this  research  was  that  of  a  motor  operating  at  rated 
speed  and  power.  This  case  was  considered  since  it  had  the  capability  of  generating  the 
largest  fault  currents.  This  model  does  not  take  into  account  how  the  magnitude  of  fault 
currents  affects  breaker  tripping  speeds.  The  effects  of  internal  faults  for  motors  operating 
at  speeds  lower  than  rated  speed  should  be  investigated. 
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The  fault  model  does  not  consider  reignition  voltage  of  the  arc.  When  the  arc  is 
extinguished  deionization  will  reduce  the  conductivity  of  the  column.  As  the  voltage 
reverses,  the  conductivity  that  existed  just  before  the  instant  the  arc  current  goes  to  zero 
must  be  reestablished.  If  the  deionization  has  been  rapid,  the  required  voltage  to 
reestablish  the  arc  will  be  consequently  higher  than  the  arc  burning  voltage.  Reignition 
voltages  can  be  of  the  order  of  0.2  to  0.8  per  unit,  Appendix  A  [2].  Introducing  the 
reignition  voltage  will  result  in  a  more  realistic  model  and  possibly  in  lower  energy 
dissipation  in  the  arc  as  the  arc  will  not  burn  as  frequently  as  in  the  current  model. 

The  model  used  does  not  account  for  the  inertia  of  the  propeller  and  the  shaft. 
This  added  inertia  will  cause  the  windmilling  shaft  to  coast  down  at  a  lower  rate.  In  this 
case  high  internal  voltages  will  be  sustained  for  a  longer  time,  thus  sustaining  the  electrical 
arc  longer. 

5.3  Suggestions  for  Future  Research 

The  flexibility  added  by  object-oriented  programming  makes  it  easier  to  add 
components  to  the  model  such  as  induction  motors  and  other  ship's  electrical  loads.  This 
could  lead  to  a  simulation  of  the  entire  ship's  electrical  distribution  system. 

For  electric  drive  ships,  the  propulsion  motors  can  be  modeled  and  incorporated 
into  the  simulation  model. 

The  fault  model  could  be  improved  to  include  the  effects  of  arc  reignition  voltages. 
The  arc  reignition  voltage  is  several  times  greater  than  the  burning  voltage  of  the  arc. 
Inclusion  of  these  voltages  into  the  model  will  provide  a  more  realistic  model  of  the  arc 
and  a  better  estimate  of  the  power  dissipated  in  the  arc. 
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The  damage  caused  by  an  arc  will  occur  in  a  very  short  period  of  time.  The  use  of 
sensors  inside  the  motor  that  can  detect  the  formation  of  a  plasma  column  should  be 
considered    This  system  would  disconnect  the  motor  immediately  upon  the  formation  of 
the  arc 

5.4  Conclusions 

This  research  has  shown  that  internal  faults  such  as  electrical  arcs  can  generate 
large  currents  and  dissipate  large  amounts  of  power  inside  a  motor.  Because  the  field  flux 
is  constant  in  permanent  magnet  machines,  this  type  of  machine  will  continue  to  generate 
an  internal  voltage  while  it  is  windmilling.  This  internal  voltage  can  continue  to  support 
the  fault  process  and  could  cause  further  damage  to  the  machine  after  it  is  disconnected 
from  its  electrical  power  source. 

For  a  large  machine  such  as  the  one  used  for  this  research,  average  powers  of 
approximately  two  to  three  megawatts  were  dissipated  at  the  arc.  Such  large  power 
dissipation  in  a  localized  spot  can  be  expected  to  cause  extensive  damage  to  the  motor. 

Shorting  the  windings  of  the  motor  was  shown  to  reduce  the  amount  of  power 
dissipation  at  the  arc.  However,  after  the  windings  are  shorted  a  significant  amount  of 
power  will  continue  to  be  dissipated  at  the  arc.  To  stop  the  process  the  machine  needs  to 
be  slowed  down  to  where  the  internal  voltage  generated  can  no  longer  support  the  arcing 
process.  Preferably  the  machine  should  be  stopped. 

Permanent  magnet  motors  are  being  considered  for  ship  propulsion  in  electric  drive 
ships.  This  research  has  shown  that  an  internal  fault  can  cause  significant  damage  to  a 
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windmilling  permanent  magnet  motor.  Further  investigation  in  the  area  of  protection 
systems  for  propulsion  permanent  magnet  motor  should  be  considered. 
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Appendix  A.  Motor  Design  Spreadsheet 


Input  Parameters 


No.  of  phases 

L/D 

Pout 

Rotor  speed,  rpm 

Frequency.  Hz 

Air  gap  size,  m 

Terminal  Voltage.  V 

Power  factor 

ks 

Ja .  Amps/mA2 

Shear,  psi 

Bm.  T  rms 

Bg.  T  rms 

Byoke,  T  rms 

Bt.  T  rms 

Magnet  Fraction 

Space  Factor 

q.  slots/pole-phase 

Conductors/slot 


3 

0.229 

40,000 

200 

60 

0.005 

1000 

0.8 

0.9 

4.00E+06 

7.5 

0.4 

0.6 

1.2 

1.08 

0.5 

0.5 

1 

1 


0) 

377 

r/s 

p 

18 

pole  pairs 

Torque 

1  42E+06 

Nt-m 

K 

9.57E+04 

A/m 

R 

2.12 

m 

B, 

0.54 

T 

Y 

0.64 

radians 

Is 

15060.0587 

Amps 

|  3  VI| 

45.18 

MVA 

Stator  Parameters 

A.  p 

0.5 

Slot  height 

5.32 

cm 

Ns 

15.71 

turns/phase 

Slot  width 

6.18 

cm 

N  slots 

108.00 

Tooth  width 

6.18 

cm 

Back  iron 

5.91 

cm 

Eff  cond  length 

1.44 

m 

Magnet  Size  Calculations 


w 


e 

Y 
s 


Reactances,  per  unitized  to  Eaf 


37.074 

cm 

0.087 

radians 

Xad 

0.3953 

p.u. 

18.531 

cm 

X 

aq 

2.3128 

p.u. 

0.087 

radians 

Zbase 

0.0335 

ohms 

13.903 

cm 

Ra 

0.0005 

ohms/phase 

0.071 

Ta 

0.0145 

p.u. 

0.16 

cm 

H 

0.1227 

sec 

Constants,  units  as  specified 

Conductivity  of  Copper.  S/m 
Density  of  copper.  Kg/mA3 
Density  of  Steel,  Kg/mA3 
Permeability  of  free  space.  H/m 
Density  of  Aluminum.  Kg/mA3 


Stator  Mass 

5.70E+07 

Teeth 

2723.35 

Kg 

8.95E+03 

Back  iron 

6227.63 

Kg 

7.80E+03 

Copper 

2293.83 

Kg 

1.26E-06 

Stator  Mass 

11244.81 

Kg 

2.70E+03 

Stator  Weight 

11 

lton 
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Appendix  A.  Motor  Design  Spreadsheet 


Calculation  of  v  (  v=  terminal  voltage/Eaf) 


v 

2.137 
2.143 
2.140 
2.142 


1.721 
1.727 
1.724 
1.726 


Magnetizing 

0.649 

Slot  leakage 

0.01 

larmonic  Leakage 

0.016 

Leakage 

0.03 

Laao 

0.433 

Xa(j 

yv  aq 

Synchronous 

0.675 

xd 

Xq 

Losses 

Mwatts 

Copper  Losses 

0.33 

Mechanical  Losses 

0.000 

Tooth  Losses 

0.00 

Core  losses 

0.011 

No-load  losses 

0.01 

Load  losses 

0.003 

Total  losses 

0.36 

n 

98.85% 

Rotor  Mass 

Magnets 

7028.819 

Poles 

6573.10 

Flux  Barriers 

10994.95 

Inner  Structure 

3784.432 

Rotor  Weight 

27.938 

5 

id 

eaf 

vnew 

0.301 

0.379 

0.994 

2.143 

0.300 

0.378 

1.003 

2.140 

0.301 

0.378 

0.999 

2.142 

0.300 

0.378 

1.001 

2.141 

ance.  ohms 

Reactance,  pu 

0.245 

3.69 

0.004 

0.05 

0.006 

0.09 

0.010 

0.15 

0.163 

2.46 

0.012 

0.366 

0.072 

2.140 

0.254 

7.594 

0.022 

0.511 

2.357 

2.285 

pu. 

0.01 

0.000 
0.0001 
0.0002 
0.0003 
0.0001 

0.01 


Kg 
Kg 
Kg 
Kg 
lton 
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Appendix  B.  Node  and  Network  Pre-Processor 


This  program  is  used  to  format  the  two  input  files,  *.in,  that  contain  the  node  and 
line  information  for  the  simulation  program.  It  is  the  first  program  that  needs  to  be 
executed.  The  output  of  this  program  are  two  files,  *  lfi  which  are  used  by  the  load  flow 
calculation  program,  Appendix  C. 

The  original  pre-processor  program  in  this  appendix  was  written  by  Professor 
James  L.  Kirtley    The  code  was  modified  by  F  R.Colberg  allow  connection  of  a 
synchronous  motor  to  a  network 

#include  <stream.h> 
#include  <stdlib.h> 

main(int  argc,  char  *argv[]) 

{ 

void  concat(char*,  char*,  char*), 
char  node_in_name[14]; 
char  line_in_name[14]; 
char  node_out_name[14]; 
char  line_out_name[14]; 

concat(argv[l],  "in",  nodeinname), 
concat(argv[2],  "in",  line_in_name); 
concat(argv[l],  "lfi",  nodeoutname); 
concat(argv[2],  "lfi",  line_out_name); 

filebuf  fD,  fl,  £2,  f3; 

if  (f0.open(node_in_name,  input)  =  0)  //  Bus  data  here 

{cout  «  "Can't  Open  "  «  nodeinname  «  "!\n"; 
exit(0); 

} 

if  (fl  open(line_in_name,  input)  =  0)  //  Line  data  here 

{cout  «  "Can't  Open  "  «  lineinname  «  "!\n"; 
exit(0); 
} 

f2.open(node_out_name,  output);  //  For  processed  node  info 

f3  open(line_out_name,  output);  //  Put  processed  line  info 

istream  nodein  (&f0), 
istream  linein  (&fl); 
ostream  nodeout  (&f2); 
ostream  lineout  (&f3); 
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int  nlines,  nnodes; 
nodein  »  nnodes; 
linein  »  nlines; 

int  nl [nlines],  n2[nlines],  line_count=0,  node_count=0; 
double  r[nlines],  xl  [nlines],  xO[nlines], 
char  sfnlines]; 


for  (int  i=0;  i<nlines,  i++) 

{ 

linein  »  nl[i]; 

linein  »  n2[i]; 

linein  »  r[i], 

linein  »  xl[i], 

linein  »  x0[i]; 

linein  »  s[i], 

if  (s[i]  ==  'c')  line_count++; 
} 


//  Get  Basic  Line  Info 


// 


This  many  to  output  file 


char  c; 

char  t[nnodes]; 

char  gname[nnodes][10]; 

int  nodefnnodes]; 

double  datl [nnodes]; 

double  dat2[nnodes]; 

int  gencount  =  0; 

for  (int  j=0;  j<nnodes;  j++) 

{ 
nodein  »  c; 

tUl  =  c; 

switch(c) 

{ 

case  'n': 

node_count++, 
break; 
case  'i': 

node_count++; 

nodein  »  dat  1  [j]; 

nodein  »  dat2[j], 

break; 
case  'd': 

break; 


//  Network  Node 


//  Voltage  Source  Node 


//  used  for  datum  Node 
//  not  in  node  count 
//  generator  node 
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case  'p': 

node_count++; 

datl[j]  =  dat2[j]=0; 

break; 
case  'g':  //  generator  or  motor  data 

nodein  »  gname[j]; 

nodein  »  node[j]; 

nodein  »  dat  1  [j  ] ; 

nodein  »  dat2[j]; 

gen_count++, 

break, 
default:  //  unrecognizable  code 

cout  «  "INPUT  FILE  ERROR  UN  RECOGNIZED  NODE  CODE'  "; 

cout«t[j]«"\n"; 

exit(O), 

} 
} 
f0.close(); 
fl.close(); 

//  First  output  file  is  line  data  for  load  flow 

lineout  «  line_count  «  "\n";  //  Top  Key  to  line  file 

for  (i=0;  i<line_count,  i++) 
if  (s[i]  =  V)  lineout  «  nl[i]  «  "  "  «  n2[i]  «  "  " 
«r[i]«"  "«xl[i]«H\n"; 

//  Second  output  file  is  node  data  for  load  flow 
//  First,  must  fix  up  generator  buses  for  actual  load 

for  (j=0;  j<nnodes,  j++) 

{ 
if(t[j]  ==  'g')  //  Generator 

switch(t[node[j]])  //  Check  target  node 

{ 

case  V:  //  Allowed  cases 

dat  1  [node[j]]  +=  dat  1  [j],       //  set  voltage 
if(dat2[node[j]]  =  0) 
dat2[node[j]]  =  dat2[j], 
break; 


i„i. 


case  p 

datl[node[j]]  +=  datlQj];       //  P  and  Q  add 
dat2[node[j]]  +=  dat2[j], 
break; 
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default: 
case  'n': 
break; 

cout  «  "INPUT  FILE  ERROR  MACHINE  ASSIGNED  TO  WRONG 
NODE!\n", 

exit(O); 

} 
} 

//  Now  we  can  actually  write  the  output  file 

nodeout  «  nodecount  «  "\n", 
for  (j=0;  j<nnodes,  j++) 

{ 
switch  (t[j]) 

{ 

case  'n':  //  Network  node:  special  p,q 

nodeout  «  "8  6  1  0  n\n",         //  p=.8,q=.6,  v=l,  d=0 

break, 
case  'i':  //  Voltage  source 

nodeout  «  "0  0  "; 

nodeout  «  dat  1  [j]  «  "  ";        //  voltage 

nodeout  «  dat2[j]  «  "  i\n";      //  phase  angle 

break, 
case  V:  //  Generator,  fixed  voltage 

nodeout  «  dat  1  [j],  //  total  power  at  note 

nodeout  «  "  0  ";  //  this  number  will  be  ignored 

nodeout  «  dat2[j];  //  fixed  voltage  at  node 

nodeout  «  "  0  v\n",  //  ignored  and  key 

break; 
case  'p':  //  Generator,  fixed  p,  q 

nodeout  «  dat  l[j]«"  ";        // P 

nodeout  «  dat2[j],  //  Q 

nodeout  «  "  1  0  p\n";  //  V  to  start  loadflow  and  key 

break, 
default:  //  Only  'g'  or  'd'  are  left 

break; 

} 
} 
f2  close(); 
f3.close(); 
} 

void  concat  (char  *  string  1,  char*string2,  char*  string) 
{ 
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for  (int  j=0;  j<strlen( string  1 );  j++) 

string[j]  =  string  l[j]; 
for  (int  i=0;  i<strlen(string2);  i++) 

string[i+j]  =  string2[i], 
} 
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Appendix  C.  Load  Flow  Program 

The  files  *  lfi  generated  by  the  pre-processor  program  are  used  by  this  program  to 
calculate  node  voltages,  line  currents  and  power  flow    This  program  generates  two  output 
file  *.lfo    These  output  files  contain  the  load  flow  information  for  the  connected  system 
The  actual  load  flow  calculations  are  performed  by  the  program  loadflow.h.  Appendix  C-l 

This  program  was  written  by  Professor  James  L  Kirtley  The  program  was 
modified  by  F.R.Colberg  to  perform  load  flow  calculations  for  nodes  and  lines  with 
motors  connected  to  them 

#include  <stream.h> 
include  "loadflow  h" 

main(int  argc,  char  *argv[] 

{ 

void  concat(char*,  char*,  char*); 
char  node_in_name[14]; 
char  line_in_name[14]; 
char  node_out_name[14], 
char  line_out_name[14]; 

concat(argv[l],  "lfi",  nodeinname); 
concat(argv[2],  "lfi",  lineinname); 
concat(argv[l],  "lfo",  nodeoutname); 
concat(argv[2],  "lfo",  lineoutname), 

filebuffO,  fl,f2,  f3, 

if  (fD.open(node_in_name,  input)  =  0)  //  Bus  data 

{cout  «  "Can't  Open  "  «  nodeinname  «  "!\n"; 
exit(0); 
} 

if  (fl  open(line_in_name,  input)  ==  0)  //  Line  data 

{cout  «  "Can't  Open  "  «  lineinname  «  "!\n"; 
exit(0); 
} 

f2.open(node_out_name,  output);  //  Processed  node  information 

f3.open(line_out_name,  output);  //  Processed  line  information 

istream  nodein  (&f0), 
istream  linein  (&fl), 
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ostream  nodeout  (&f2); 
ostream  lineout  (&f3), 

const  double  crit=le-7; 

cnetn(fl,  fD); 

fD.closeO; 

fl.close(); 

double  e=l, 

for  (int  k=0;  k<  10000;  k++) 

if((e  =  n  improvevoltagesO)  <  crit)  break, 
cout  «  "That  Took  "  «  k  «  "  Iterations\n"; 

n.report_node_voltages(f2); 
n  report_line_currents(f3 ) , 

f2.close(); 
f3  closeQ; 


} 


void  concat  (char*  string  1,  char*string2,  char*  string) 

{ 

for  (int  j=0;  j<strlen(stringl);  j++) 

string[j]  =  stringl|j]; 

for  (int  i=0,  i<strlen(string2),  i++) 

string[i+j]  =  string2[i], 

} 
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Appendix  C-l.  Load  flow  calculation  program 

This  program  solves  the  load  flow  problem  for  the  system  defined  by  the  files  *  in. 
This  code  is  part  of  Appendix  C.  The  original  program  written  by  Professor  James  L 
Kirtley  was  modified  by  F  R  Colberg  to  add  destructor  functions  and  to  delete  reference 
to  voltage  regulator  files 

include  "clinch" 
#include  "cnode.h" 
#include  <stream  h> 
^include  <libc  h> 
#include  <stdlib  h> 
#include  "Complex  h" 
#define  Complex  complex 

class  cnet  { 
int  nlines;  //  number  of  lines 

int  nnodes;  //  number  of  nodes  in  network 

cline  *lptr,  //  pointer  to  lines 

cnode  *nptr;  //  pointer  to  nodes  (buses) 

int  *node_incidence; 
Complex  *node_admittance; 

public: 
cnet  (filebuf  ft),  filebuf  fl )  //  Build  network  from  files  *  lfi 

{ 

int  nl,  n2; 
char  t, 
double  r,  x; 
double  p,  q,  v,  d; 

istream  from  (&f0),  //  This  file  contains  line  information 

from  »  nlines, 

lptr  =  new  cline[nlines];  //  Build  the  lines  defined  in  input  file 

for  (int  i=0;  i<nlines;  i++) 

{ 

from  »  nl ,  //  Entry  node 

from  »  n2,  //  Exit  node 

from  »  r;  //  Line  resistance 

from  »  x;  //  Line  reactance 

lptr[i].set_cline(nl,  n2,  r,  x); 
} 

for  (i=0;  i<nlines,  i++) 
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lptr[i]  rept(i); 

istream  fh  (&fl );  //  This  file  contains  the  node  information 

fh  »  nnodes, 

nptr  =  new  cnodejnnodes],    //  Build  the  nodes 

for  (int  j=0;  j<nnodes;  j++) 

{ 

fn  »  p;  //  Initial  node  real  power 

fn  »  q;  //  Initial  node  reactive  power 

fh  »  v,  //  Initial  node  voltage  magnitude 

fh  »  d,  //  Initial  node  angle 

fn  »  t,  //  Node  type 
nptr[j].set_cnode(p,  q,  v,  d,  t), 
} 

for  (j=0;  j<nnodes;  j++) 

nptrfj]  rept(j); 
nodeincidence  =  new  int[nnodes*nlines]; 

for  (j=0;  j<nnodes;  j++) 
for  (i=0;  i<nlines;  i++) 
node_incidence[i+j*nlines]  =  0, 

for  (i=0;  i<nlines;  i++)       //  Build  node-incidence  matrix 

{ 

node_incidence[lptr[i]  report_node_a()*nlines+i]  =  1; 
node_incidence[lptr[i].report_node_b()*nlines+i]  =  -1; 

} 

//  Build  the  node-admittance  matrix 

nodeadmittance  =  new  Complex[nnodes*nnodes]; 

for  (j=0;  j<nnodes,  j++)  //  Allocate  space  for  admittance  matrix 

for  (i=0;  Knnodes;  i++) 
node_admittance[j*nnodes+i]  =  (0,0), 

Complex  *tmat  =  new  Complex[nnodes*nlines], 

for  (i=0;  i<nlines;  i++) 
for  (j=0,  j<nnodes;  j++) 
tmat[i*nnodes+j]  =  lptr[i].admit()*node_incidence[j*nlines+i]; 

for  (j=0;  j<nnodes;  j++) 
for  (i=0;  i<nnodes,  i++) 
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for  (int  k  =  0;  k<nlines;  k++) 
node_admittance[j*nnodes+i]  += 
node_incidence[j*nlines+k]*tmat[k*nnodes+i], 

delete  [nnodes*nlines]tmat, 

}    //  End  of  network  construction  step  cnet() 

Complex  node_current(int  i)  //  Current  at  node  i 

{ 

Complex  I; 
I  =  (0,0), 
for  (int  j  =  0,  j<nnodes,  j++)  { 

I  +=  nptr[)].node_voltage()  *  node_admittance[i*nnodes+j];} 
return(I); 

} 

void  report_node_power() 

x 
Complex  I,  P; 

for  (int  j  =  0;  j<nnodes;  j++) 

{ 

I  =  (0,0); 

for  (int  i=0;  i<nnodes;  i++) 

I  +=  nptr[i].node_voltage()  *  node_admittance[j*nnodes+i]; 
P  =  nptr[j].node_voltage()  *  conj(I); 
cout  «  "Node  "  « j  «  "  P  +  j  Q  =  "  «  P  «  "\n", 
} 
} 

double  improve_voltages()  //  One  step  in  voltage  improvement 

{ 

Complex  C,  Y, 

double  e  =  0,  //  error  accumulator 

for  (int  i=0;  Knnodes;  i++) 

{ 

C  =  nodecurrent(i); 

Y  =  node_admittance[i+nnodes*i]; 

e  +=  nptr[i]  improve_voltage(C,  Y)/nnodes, 

} 
return(e); 

} 

void  report_node_voltages(filebuf  £2)  //Output  node  voltages 

{ 
Complex  V,  I,  P, 
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ostream  to(&f!2); 

for  (int  i=0;  i<nnodes,  i++) 

{ 

V  =  nptr[i].node_voltage(), 

I  =  nodecurrent(i); 

P  =  V  *  conj(I); 

to  «  abs(V)  «  "  "  «  arg(V)  «  "  " 

«  real(P)  «  "  "  «  imag(P)  «  "V; 
} 
} 
void  report_line_currents(filebuf  f3)  //Output  line  currents 

{ 

Complex  c,  vl,  v2,  I; 
int  nl,  n2; 
ostream  to(&f3); 
for  (int  i=0;  Knlines,  i++) 

{ 

c  =  lptr[i].admit(); 
nl  =  lptr[i].report_node_a(); 
n2  =  lptr[i]  report_node_b(); 
vl  =  nptr[nl].node_voltage(); 
v2  =  nptr[n2].node_voltage(), 
I  =  c  *  (vl  -  v2); 
to  «  abs(I)  «  "  "  «  arg(I)  «  "\n", 

} 
} 

~cnet();  //Destructor  function 

} ;  //  End  of  definition  of  class  cnet 

cnet::~cnet()  { 

delete  [nlines]lptr; 

delete  [nnodes]nptr, 

delete  [nnodes*nnodes]node_admittance, 

delete  [nnodes*nlines]node_incidence, 
} 


Appendix  C-2.  Node  Voltage  Calculation 

This  file  calculates  node  voltages  for  the  load  flow  program    It  is  part  of  Appendix 
C  and  needs  to  be  compiled  as  part  of  it.  This  file  was  written  by  Professor  James  L. 
Kirtley 

include  "Complex.h" 
#include  <stream.h> 
#include  <stdlib  h> 
#define  complex  Complex 

class  cnode  { 

double  p,  q,  v,  d;  //  Real,  reactive  power,  voltage,  angle 

char  t,  //  Node  type 

public: 

void  setcnode  (double  pwr,  double  qwr,  double  vt,  double  delt,  char  type) 

{ 

p  =  pwr; 

q  =  qwr, 

v  =  vt; 

d  =  delt, 

t  =  type; 
} 

Complex  node_voltage()        //  Report  voltage  at  node  i 

{ 
return(polar(v,  d)), 

} 

double  improve_voltage(Complex  I,  Complex  Y)  //  One  step  in  a  Gauss-Seidel 
{  // 1  is  current  from  network 

Complex  cv,  vo, 
double  qt,  e,  a=1.6; 
switch(t) 

{ 

case  V:  //  Infinite  bus 

e  =  0; 

break; 
case  'p':  //  Defined  p-q  node 

vo  =  polar(v,  d); 

cv  =  -Complex(p,-q)/(conj(vo)*Y)+  I/Y,  //  Voltage  correction 

vo  +=  a*cv, 

v  =  abs(vo), 

d  =  arg(vo), 
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e  =  abs(Complex(p,q)-vo*conj(I)), 

break; 
case  V:  //  Defined  p  and  |v| 

vo  =  polar(v,d); 

qt  =  imag(vo*conj(I)),  //  Calculated  value  of  q 

cv  =  Complex(p,-qt)/(conj(vo)*Y)-  I/Y, 

vo  +=  a*cv;  //  First-order  correction  to  vo 

vo  *=  v/abs(vo), 

d  =  arg(vo), 

e  =  abs(p  -  real(vo*conj(I))), 

cout  «  "vo  =  "  «  vo  «  "  I  =  "  « I «  "  e  =  "  «  e  «  "\n' 

break; 
default: 

cout  «  "Incorrect  Node  Type  Used  in  Input  \n"; 

exit(O), 

} 
return(e); 

} 

void  rept(int  i) 

{ 

cout  «  "  Node  "  «  i  «  "  Type  "  «  t  «  "  P=  "  «  p  «  "  Q=  " 
«  q  «  "  V=  "  «  v  «  "  Delt=  "  «  d  «  "\n", 

} 
Complex  report_target_power() 

{ 
return(Complex(p,q)), 

} 

}; 

//  End  of  definition  of  class  node 


83 


Appendix  C-3.  Line  Admittance  Calculation 

This  file  calculates  the  line  admittances    It  is  part  of  Appendix  C  the  load  flow 
program.  This  file  was  written  by  Professor  James  L  Kirtley 

#include  "Complex. h" 

typedef  class  complex  Complex; 

class  cline  { 

int  na,  nb, 

double  r,x, 
public: 

void  set_cline  (int  nodea,  int  nodeb,  double  res,  double  react) 

{ 

na  =  nodea; 

nb  =  nodeb, 
r  =  res; 
x  =  react; 

} 

Complex  admit()     //  Report  line  admittance 

{ 

Complex  z,  y; 
z  =  Complex(r,  x); 
y=  1.0/z; 
return(y); 

} 
int  report_node_a() 

{ 
return(na); 

} 
int  report_node_b() 

{ 
return(nb); 

} 
void  rept(int  i) 

{ 

cout  «  "  Line  "  «  i  «  "  from  "  «  na  «  "  to  " 
«  nb  «  "  z  =  "  «  r  «  "  +  j"  «  x  «  "\n"; 

} 

}; 

//  End  of  definition  of  class  cline 


84 


Appendix  D.  Line  Simulation  Input  File 

This  program  uses  the  power  flow  information  calculated  in  Appendix  C  to  build 
the  lines  for  the  simulation  program.  This  program  will  generate  an  output  file  *.nct 
which  contains  the  line  information  for  the  simulation  program 

Original  program  written  by  Professor  James  L.  Kirtley    The  code  was  modified 
by  F.R.Colberg  to  delete  reference  to  voltage  regulators  and  to  allow  connection  of  a 
motor  to  the  network 

#include  <stream.h> 
#include  <stdlib.h> 
#include  <math.h> 
#include  "Complex. h" 
#defme  Complex  complex 

main(int  argc,  char  *argv[]) 

{ 
const  double  al  =  2.094395 1 ;  112  pi/3 

void  concat(char*,  char*,  char*); 

char  node_in_name[14], 

char  line_in_name[14]; 

char  node_V_name[14]; 

char  line_I_name[14]; 

char  line_out_name[14]; 

char  node_output_file_name[  1 4] , 

char  line_output_file_name[14]; 

char  mot_output_file_name[14]; 

concat(argv[l],  "in",  nodein  name), 
concat(argv[2],  "in",  lineinname); 
concat(argv[l],  "lfo",  nodeVname), 
concat(argv[2],  "lfo",  linelname); 
concat(argv[2],  "net",  lineoutname), 
concat(argv[l],  "  gp",  motoutputfilename); 

filebuf  fO,  fl,  f2,  f3,  f4,  fg; 

if  (f0.open(node_in_name,  input)  =  0)  //  Bus  data 

{cout  «  "Can't  Open  "  «  nodeinname  «  "!\n", 
exit(0); 
} 
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if  (fl  open(line_in_name,  input)  =  0)  //  Line  data 

{cout  «  "Can't  Open  "  «  lineinname  «  "!\nH; 
exit(0), 

} 

if  (f2.open(node_V_name,  input)  ==  0)  //  Bus  data 

{cout  «  "Can't  Open  "  «  nodeVname  «  "!\n", 
exit(0); 

} 

if  (B  open(line_I_name,  input)  ==  0)  //  Line  data 

{cout  «  "Can't  Open  "  «  linelname  «  "!\n", 
exit(0); 

} 

f4.open(line_out_name,  output),  //  Processed  line  information 

fg.open(mot_output_file_name,  output); 

istream  nodein  (&f0), 
istream  linein  (&fl); 
istream  nvin  (&f2); 
istream  liin  (&f3); 
ostream  lineout  (&f4); 
ostream  motout  (&fg); 

int  nlines, 
linein  »  nlines, 

int  nl  [nlines],  n2 [nlines]; 

double  rfnlines],  xl [nlines],  x0[nlines], 

char  sfnlines]; 

int  lcode[nlines]; 

for  (int  i=0,  Knlines;  i++)  //  Get  Basic  Line  Info 

{ 

linein  »  nl [i], 

linein  »  n2[i]; 

linein  »  r[i]; 

linein  »  xl[i], 

linein  »  x0[i], 

linein  »  s[i], 
} 

linein  » lineoutputfilename; 
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int  nnodes; 
nodein  »  nnodes; 

char  c; 

char  t[nnodes]; 

char  mname[nnodes][10]; 

int  node[nnodes], 

double  datl [nnodes]; 

double  dat2[nnodes]; 

int  gencount  =  0, 

double  nvm[nnodes],  nva[nnodes],  p[nnodes],  q[nnodes], 

cout  «  "Reading  Network  Data\n"; 

for  (int  j=0,  j<nnodes;  j++) 

{ 
nodein  »  c, 

tOJ  =  c, 

cout  «  "j  =  "  « j  «  "  c  =  "  «  c  «  "\n", 

switch(c) 

{ 
case  'n': 

nvin»nvm[j];  //  Network  Node 

nvin  »  nva[j]; 

nvin  »  p[j], 

nvin  »  q[j]; 

cout  « j  «  "  n  nvm  ="  «  nvm[j]  «  "  nva  ="  «  nva[j] 

«  »  p  =••  «  p|j]  «  -  q  ="  «  qrj]  «  -in- 
break; 
case  V:  //  Voltage  Source  Node 

nodein  »  dat  1  Qj], 
nodein  »  dat2[j]; 
nvin  »  nvm[j]; 
nvin  »  nva[j]; 
nvin  »  p[j]; 
nvin  »  q[j], 

cout  «j  «  "  i  nvm  ="  «  nvm[j]  «  "  nva  ="  «  nvajj] 
«  »  p  =••  «  prj]  «  •■  q  =••  «  q[j]  «  "v; 

break; 
case  'd':  //  used  for  datum  Node 

break,  //  not  in  node  count 

case  V:  //  generator  node 

datl[j]  =  dat2[j]=0; 

nvin  »  nvm[j]; 

nvin  »  nva[j]; 
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} 


nvin  »  p[j]; 

nvin  »  q[j]; 

cout  « j  «  "  v  nvm  ="  «  nvm[j]  «  "  nva  ="  «  nva[j] 
«  "  p  ="  «  p[j]  «  "  q  ="  «  q[j]  «  »\n"; 

break; 
case  'p': 

datl[j]  =  dat2[j]=0; 

nvin  »  nvm[j]; 

nvin  »  nva[j], 

nvin  »  p[j]; 

nvin  »  q[j]; 

cout  « j  «  "  p  nvm  ="  «  nvm[j]  «  "  nva  ="  «  nva[j] 
«  "  p  ="  «  p[j]  «  "  q  ="  «  q[j]  «  "\n", 

break; 
case  'g': 

nodein  »  mname[j]; 

nodein»  node[j]; 

nodein  »  datl[j]; 

nodein  »  dat2[j]; 

gen_count++; 

break; 
default:  //  unrecognizable  code 

cout  «  "INPUT  FILE  ERROR:  UNRECOGNIZED  NODE  CODE!  ", 

cout  «  t[j]  «  "\n"; 

exit(O); 

break; 
} 


//  motor  data 


unrecognizable  code 


cout «  "Node  Data\n"; 
for  (j=0;  j<nnodes;  j++) 

cout  «  "Node  "«j«"\tv=  "«nvm[j]«"\ta=  "«nva[j]«"\n'\ 

int  ntlines  =  nlines  +  gencount;  //  total  number  of  lines 

for  (i=0;  i<nlines;  i++) 
linein  »  lcode[i], 


lineout  «  ntlines  «  "\n"; 


//  start  output 


double  ia,  ib,  ic; 
int  sa,  sb,  sc; 
double  lim,  lia, 


for  (i=0;  Knlines;  i++) 


//  ordinary  lines  first 
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{ 

if(s[i]=='o')  //open  line 

{ 

ia  =  ib  =  ic  =  0, 
sa  =  sb  =  sc  =  0, 

} 
else 

I 
liin  »  lim, 

liin  »  lia, 

ia  =  lim  *  cos  (lia); 

ib  =  lim  *  cos  (lia-  al); 

ic  =  lim  *  cos  (lia  +  al), 

sa  =  sb  =  sc  =  1 ; 

} 
lineout  «  nl[i]  «  "  "  «  n2[i]  «  "  " 

«  r[i]  «  "  "  «  xl[i]  «  "  "  «  x0[i]  «  "  " 
«  sa  «  "  "  «  sb  «  "  "  «  sc  «  "  " 
«  ia  «  "  "  «  ib  «  "  "  «  ic  «  "  " 
«  "\n"; 
} 

fD.close(); 
fl.closeO; 

f2.close(), 
f3.close(); 

//Build  the  machine  connection  stubs 

double  gxl,  gxO,  ra,  tdopp,  tqopp, 
double  xd,  xq,  xdpp,  xqpp,  ta,  h,  omz,  qfract; 
double  P,  Q; 
complex  I,  V; 

filebuffS; 

i=0, 
int  k; 

for  (j=0;  j<nnodes,  j++) 
{ 

if(tD]=='g') 
{ 
if  (f5. open  (mname[j],  input)  ==  0) 

{ 
cout  «  "Missing  Machine  File  "  «  mname[j]  «  "\n"; 
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exit(O); 

} 
istream  gin  (&f5), 

gin  »  xd; 

gin  »  xq; 
gin  »  xdpp, 
gin  »  xqpp, 
gin  »  gxO, 
gin  » tdopp, 
gin  » tqopp; 
gin  » ta; 
gin  »  h; 
gin  »  omz, 


f5.close(); 


gxl  -  0.5  *  (xdpp  +  xqpp),  //  positive  sequence  reactance 

ra  =  gxl  /  (omz  *  ta),  //  line  resistance 

//  Have  to  get  current  for  the  machine 


switch(t[node[j]]) 

{ 

case  'p': 
P  =  datl[j], 

Q  =  dat2[j]; 
break; 
case  V: 

P  =  datl[j], 

qfract=  1 ; 

for  (k=0;  k<nnodes,  k++) 

( 

if  ((k<j)  &&  (t[k]=='g')  &&  (node[k]  ==  node[j])) 

{ 

qfract  =  dat2[j], 
break, 

} 
else  if  ((k>j)  &&  (t[k]  =  'g')  &&  (node[k]  =  nodefj])) 
qfract  -=  dat2[k]; 

} 
Q  =  qfract*q[node[j]], 
break; 
default 


90 


cout  «  "  Apparent  Generator  At  Non-Generator  Node!\n"; 
exit(O); 

} 
I  =  conj(complex(P,Q)/polar(nvm[node[j]],  nva[node[j]])), 

cout  «  "Generator  Stub:  P  =  "  «  P  «  "  Q  =  "  «  Q  «  H\n"; 
cout  «  "Node  "  « j  «  "  To  "  «  node[j]  «  "  I  =  "  «  I  «  "\n", 
cout  «  "Node  v  =  "  «  nvm[node[j]]  «  "  Angle  =  "  «  nva[nodeO]] 

«  "\n"; 
motout  «  nvm[node[j]]  «  "  "  «  nva[node[j]]  «  "  " 

«  p  «  "  "  «  Q  «  "  "  «  gxl  «  "  "  «  ra  «  "\n", 

lim  =  abs(I); 
lia  =  arg(I), 

cout  «  "I  =  "  «  lim  «  "\t Angle  =  "  «  lia  «  "\tra  =  " 
«  ra  «  "\tgxl  =  "  «  gxl  «  "\n"; 

ia  =  lim  *  cos  (lia); 
ib  =  lim  *  cos  (lia  -  al); 
ic  =  lim  *  cos  (lia  +  al); 
sa  =  sb  =  sc  =  1 , 

lineout  «j  «  "  "  «  node[j]  «  "  " 
«  ra  «  "  "  «  gxl  «  "  "  «  gxO  «  "  " 
«  sa  «  "  "  «  sb  «  "  "  «  sc  « 
«  ia  «  "  "  «  ib  «  "  "  «  ic  « 
«  "\n"; 

} 
} 
lineout  «  lineoutputfilename  «  "\n", 
for  (i=0;  Knlines,  i++) 

lineout  «  lcode[i]  «  "\n"; 
for  (i=0;  i<gen_count,  i++) 
lineout  «  "7\n", 
f4.close(); 
fg.close(); 

} 

void  concat  (char  *stringl,  char*string2,  char*string) 

{ 
for  (int  j=0;  j<strlen(stringl),  j++) 

string£j]  =  string  l[j]; 
for  (int  i=0;  i<strlen(string2),  i++) 

string[i+j]  =  string2[i], 
} 
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Appendix  E.  Synchronous  Machine  /  Network  Simulator 

This  program  is  the  actual  simulation  program    It  uses  the  output  files  generated 
by  Appendices  B  through  D  and  the  file  called  base  to  run  the  simulation    The  output 
from  the  simulation  goes  to  the  output  files  designated  in  the  *  in  and  base  files 

Program  written  by  Professor  James  L.Kirtley,  modified  by  F  R  Colberg  to 
incorporate  the  permanent  magnet  motor  model  into  the  simulation. 

#include  <stream  h> 

#include  <stdlib  h> 

#include  "Complex. h" 

include  "motl.h" 

typedef  class  complex  Complex, 

main(int  argc,  char  *argv[]) 

{ 

double  dt,  omz; 
int  np,  no; 
const  double  al  =  2.0943951;  //  2  pi/3 

void  concat(char*,  char*,  char*); 
char  node_in_name[14]; 
char  node_V_name[14]; 
char  line_in[14]; 
char  mot_in_name[14]; 

concat(argv[l],  "in",  nodeinname); 
concat(argv[l],  "  lfo",  nodeVname); 
concat(argv[l],  "gp",  motinname); 
concat(argv[2],  "net",  linein); 

filebuffO,  fl,£2,  f3,  fp, 

if  (ft)  open(node_in_name,  input)  ==  0) 
{cout  «  "Can't  Open  "  «  argv[l]  «  "!\n"; 
exit(0), 
} 

if  (fl .  open(node_V_name,  input)  =  0)  //  Bus  data  here 

{cout  «  "Can't  Open  "  «  nodeVname  «  "!\n", 
exit(0); 
} 
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if  (f2.open(line_in,  input)  =  0)  //  Line  data 

{cout  «  "Can't  Open  "  «  linein  «  "!\n"; 
exit(0); 

} 

if  (f3.open(argv[3],  input)  ==  0)  //  Base  data 

{cout  «  "Can't  Open  "  «  argv[3]  «  "!\n"; 
exit(0); 

} 

if  (fp.open(mot_in_name,  input)  =  0)  //  Machine  data 

{cout  «  "Can't  Open  "  «  motinname  «  "!\n", 
exit(0), 

} 

istream  nodein  (&f0); 
istream  nvin  (&fl), 
istream  linein  (&f2), 
istream  basein  (&f3); 
istream  motin  (&fp), 

cout  «  "\n  Getting  Simulation  Data  From  "  «  argv[3]  «  "\n" 

basein  »  omz, 
basein  »  dt, 
basein  »  np; 
basein  »  no, 

int  nevents; 
basein  »  nevents; 
double  event_time[nevents]; 
int  event_line[nevents], 
char  event_type[nevents]; 

cout  «  "  omz  =  "  «  omz  «  "  dt  =  "  «  dt  «  "  no  =  "  «  no 

«  "  np  =  "  «  np  «  "\n", 

cout  «  "  There  are  "  «  nevents  «  "Events An"; 

for  (int  ne=0,  ne<nevents,  ne++) 

{ 

basein  »  event_time[ne], 
basein  »  event_line[ne], 
basein  »  event_type[ne], 

} 
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for  (ne=0;  ne<nevents,  ne++) 
cout  «  "  At  time  =  "  «  event_time[ne]  «  "  Line  =  " 
«  event_line[ne]  «  "  Type  =  "  «  event_type[ne]  «  "\n"; 

f3.close(); 


int  nlines, 
linein  »  nlines, 

int  nl,  n2,  sa,  sb,  sc; 
double  r,  xl,  xO,  ia,  ib,  ic; 

line  *lptr; 

lptr  =  new  line[nlines]; 

cout  «  "\n  Getting  Line  Data  From  "  «  linein  «  "\n"; 
cout  «  "  There  are  "  «  nlines  «  "  LinesW; 

for  (int  i=0;  i<nlines;  i++)  //  Here  we  set  up  the  lines 

{ 

linein  »  nl, 

linein  »  n2; 

linein  »  r, 

linein  »  xl; 

linein  »  xO, 

linein  »  sa; 

linein  »  sb, 

linein  »  sc; 

linein  » ia; 

linein  »  ib; 

linein  »  ic; 

cout  «  "Line  "  «  i  «  "  Between  Nodes  "  «  nl  «  "  and  "  «  n2  «  "\n"; 

cout  «  "xl=  "  «  xl  «  "  x0=  "  «  xO  «  "  r=  "  «  r  «  "\n"; 

lptr[i]  setline(nl,  n2,  xl,  xO,  r,  omz,  sa,  sb,  sc,  ia,  ib,  ic); 

lptr[i].setup(); 
} 

char  line_out_name[14]; 
int  lineoutkey  [nlines]; 

linein  » lineoutname; 
for  (i=0;  i<nlines,  i++) 
linein  » line_out_key[i]; 
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f2.close(), 

filebuffl; 

fl.open(line_out_name,  output), 
ostream  lineout  (&fl); 

cout  «  "Line  Output  Data  Will  Go  To  "  «  lineoutname  «  "\nM, 
cout  «  "Line  Output  Keys  Are:  "; 
for  (i=0;  Knlines,  i++) 

cout  «  line_out_key[i]  «  "  "; 
cout «  "\n"; 

cout  «  "  Brief  Line  Summary\n", 
for  (i=0;  i<nlines;  i++) 
cout  «  "  From  "  «  lptr[i].report_node_a()  «  "  To  " 
«  lptr[i].report_node_b()  «  "  pointer  "  «  &lptr[i]  «  "\n"; 


cout  «  "\n  Getting  Node  Input  Data  From  "  «  nodeinname  «  "\n"; 
cout  «  "  Getting  Load  Flow  Data  From  "  «  nodeVname  «  "\n", 

int  nnodes, 
nodein  »  nnodes, 

cout  «  "  There  are  "  «  nnodes  «  "  Nodes\n", 

char  c; 

char  t[nnodes], 

char  mname[nnodes][10], 

int  node[nnodes], 

double  datl  [nnodes]; 

double  dat2[nnodes], 

int  gencount  =  0; 

int  netnodecount  =  0; 

int  vnodecount  =  0; 

double  nvm[nnodes],  nva[nnodes],  p[nnodes],  q[nnodes], 
int  ptr[nnodes], 

for  (int  j=0;  j<nnodes;  j++) 

{ 

nodein  »  c, 

cout  «  "  Node  "  «j  «  "  c=  "  «  c  «  "\n"; 
switch(c) 
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{ 


nvin  »  nvm[j];  //  Network  Node 

nvin  »  nva[j], 

nvin  »  p[j]; 

nvin  »  q[j]; 

cout  «  "Network  (type  n)  Node\n", 

cout  « j  «  "  n  nvm  ="  «  nvm[j]  «  "  nva  ="  «  nva[j] 

«  "  p  ="  «  p|J]  «  h  q  =,.  «  q[J]  «  H^H. 

tD]  =  V; 

ptrjj]  -  netnodecount; 
net_node_count++; 
break; 
case  V:  //  Voltage  Source  Node 

nodein  »  datl[j]; 
nodein  »  dat2[j], 
nvin  »  nvm[j], 
nvin  »  nva[j], 
nvin  »  p[j]; 
nvin  »  q[j]; 

cout  «  "Voltage  Source  (type  i)  Node\n", 
cout  «j  «  "  i  nvm  ="  «  nvm[j]  «  "  nva  ="  «  nva[j] 

«    "    p   ="    «   p0]    «    H    q   =M    «   qFj-|    «    H^M. 

tD]  =  V; 

ptrjj]  =  vnodecount, 
v_node_count++; 
break, 
case'd':  // Datum  Node 

t[j]  =  V; 

ptr[j]  =  vnodecount; 
v_node_count++; 

datl[j]  =  dat2[j]  =  nvm[j]  =  nva[j]  =  0; 
cout  «  "Datum  (type  d)  Node\n"; 
break,  //  Voltage  node 

case  V:  //  Generator  node 

datl[j]  =  dat2[j]=0, 
nvin  »  nvm[j], 
nvin  »  nva[j], 
nvin  »  p[j], 
nvin  »  qjj]; 


nvin  »  qyj; 
cout  «  "Generator  (type  v)  Node\n", 
cout  « j  «  "  v  nvm  ="  «  nvm[j]  «  "  nva  ="  « 
«  «  p  =»  «  prj]  «  «  q  =••  «  q[j]  «  »\n"; 

tD]  =  W; 

ptr[j]  =  netnodecount; 


nvaD] 
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net_node_count++; 
break; 
case  'p': 
datl[)]  =  dat2[j]=0; 
break, 

nvin  »  nvm[j]; 
nvin  »  nva[j]; 
nvin  »  p[j], 
nvin  »  q[j]; 

cout  «  "Generator  (type  p)  Node\n"; 
cout  «j  «  "  p  nvm  ="  «  nvm[j]  «  "  nva  ="  «  nva[j] 
«  i-  p  =•>  «  p|j]  «  -  q  =•«  «  qjj]  «  "V; 

tD]  =  w; 

ptrjj]  =  netnodecount; 
net_node_count++; 
break; 
case'g':  // Motor  data 

nodein  »  mname[j], 
nodein  »  node[j]; 
nodein  »  datljj]; 
nodein  »  dat2[j]; 

tD]  -  'g'; 

ptr[j]  =  gencount; 
cout  «  "MotorVn"; 

cout  «  "Mot  File  "  «  mnameO]  «  "\n"; 
gen_count++, 
break; 
default:  //  Unrecognizable  code 

cout  «  "INPUT  FILE  ERROR:  UN-RECOGNIZED  NODE  CODE' 
cout  «  t[j]  «  "\n", 
exit(O); 
break; 

} 


char  node_out_file_name[14], 
int  node_out_key[nnodes], 

nodein  »  nodeoutfilename, 

for  (j=0;  j<nnodes,  j++) 
nodein  »  node_out_key[)], 

f0.close(); 
fl  close(); 
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cout  «  "  Node  Voltage  Data  Will  Go  To  File  "  «  nodeoutfilename  «  M\n' 
cout  «  "  Node  Keys  Are:  "; 
for  (j=0;  j<nnodes;  j++) 

cout  «  node_out_key[j]  «  "  "; 
cout  «  "\n"; 

filebuf  fn, 

fn.open(node_out_file_name,  output); 
ostream  nodeout  (&fn); 

cout  «  "Network  Data  Read  \n"; 

motor*  motr;  //  Pointer  to  motors 

bus*  nptr,  //  Pointer  to  network  buses 

vbus*  vptr,  //  Pointer  to  voltage  buses 


motr  =  new  motor[gen_count], 
nptr  =  new  bus[net_node_count], 
vptr  =  new  vbus[v_node_count], 

cout  «  "Space  Allocated  for  Motors  and  Buses  \n", 

//  Setup  the  various  nodes 

filebuf  fg, 
int  jg=0; 
int  jn=0; 
int  jv=0; 

line**  linebuf; 

int  locallinecount,  locallineno; 

cout  «  "Ready  to  setup  nodes\n"; 

for  (j=0;  j<nnodes;  j++) 

{ 

cout  «  "Node  "  « j  «  "  Code=  "  «  t[j]  «  "\n", 

switch(t[j]) 

{ 

case  'g': 
if  (fg.open(mname(j],  input)  =  0) 

{ 

cout  «  "Can't  Open  "  «  mname[j]  «  "!  \n", 
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exit(O), 
} 


motr[jg].genset(fg); 
motr[jg].set_node(j); 


jg++; 

fg.close(); 
break; 
case  V:  //  Voltage  source  nodes 

cout  «  "Node  "  « j  «  "  Voltage  Source  V  =  " 
«  datlQj]  «  "  Angle  =  "  «  dat2jj]  «  "\n", 

vptr[jv].setvbus(datl[j],  dat2[j]); 

jv++; 
break; 
case  'n':  //  Network  nodes 

cout  «  "--Node  "  «j  «  "\n", 

local_line_count=0; 

//  Count  connected  lines 
for  (i=0;  Knlines;  i++) 
if  ((lptr[i].report_node_a()  =  j)  || 
(lptr[i].report_node_b()  ==  j)) 

{ 

local_line_count++; 

cout  «  "  Line  "  «  i  «  "\n"; 
} 

linebuf  =  new  line  *[local_line_count], 
local_line_no=0; 

for  (i=0;  Knlines;  i++)        //  Copy  the  pointers 
if  ((lptr[i].report_node_a()  ==j)  || 
(lptr[i].report_node_b()  ==  j)) 
linebuf[local_line_no++]  =  &lptr[i], 

nptr[jn].setbus(i,  locallinecount,  linebuf), 

cout  «  "Node  "  « j  «  "  Network  Node,  "  «  locallinecount 

«  "  Lines  Connected  \n"; 
for  (i=0;  i<local_line_count,  i++) 

cout  «  "  From  "  «  linebuf[i]->report_node_a() 
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«  "  To  "  «  linebufIi]->report_node_b() 
«  "  pointer  =  "  «  linebuf[i]  «  "\n"; 
jn++; 

delete  linebuf, 
break, 
default: 
break; 

} 
} 

//  Next,  set  the  bus  pointers  for  each  of  the  lines 

bbus*  busa, 
bbus*  busb; 
int  nn, 

cout  «  "  Summary  of  Bus  Pointers\n"; 
for  (j=0;  j<nnodes,  j++) 
switch(t[j]) 

{ 

case  'g': 

cout  «  "  Node  "  «j  «  "  Type  g  "  «  "  ptr  =  "  «  ptr[j] 
«  "  pointer  =  "  «  &motr[ptr[j]]  «  "\n"; 

break, 
case  V: 

cout  «  "  Node  "  « j  «  "  Type  v  "  «  "  ptr  =  "  «  ptr[j] 
«  "  pointer  =  "  «  &vptr[ptr[j]]  «  "\n"; 

break, 
case  'n': 

cout  «  "  Node  "  « j  «  "  Type  n  "  «  "  ptr  =  "  «  ptr[j] 
«  "  pointer  =  "  «  &nptr[ptr[j]]  «  "\n"; 

break, 
default: 

break, 
} 

for  (i=0;  i<nlines;  i++) 

{ 

nn  =  lptr[i]  report_node_a(), 
switch(t[nn]) 

{ 
case  'g': 

busa  =  &motr[ptr[nn]]; 

motr[ptr[nn]].set_line(&lptr[i]); 
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break: 


case  V: 


busa  =  &vptr[ptr[nn]], 
break; 
case  'n': 
busa  =  &nptr[ptr[nn]], 
break; 

} 
nn  =  lptr[i].report_node_b(); 
switch(t[nn]) 

{ 
case  'g': 

busb  =  &motr[ptr[nn]]; 

break, 
case  V: 

busb  =  &vptr[ptr[nn]], 

break; 


case  'n': 


busb  =  &nptr[ptr[nn]]; 
break; 
} 

lptr[i] .  set_bus_pointers(busa,  busb), 
} 

//  Set  the  initial  state  values  for  the  machine(s) 

double  VM,  VA,  P,  Q,  X,  R; 

double  tm[gen_count],  eafjgencount]; 

for  (j=0,  j<gen_count;  j++) 

{ 
motin  »  VM; 
motin  »  VA; 
motin  »  P, 
motin  »  Q, 
motin  »  X; 
motin  »  R; 

cout  «  "\n  Setting  Motor  "  « j  «  "  Initial  ConditionsV, 
cout  «  "  vm=  "  «  VM  «  "  va=  "  «  VA  «  "  P=  "  «  P 
«  »  q=  h  «  q  «  .-  x=  ■■  «  x  «  "  R=  "  «  R  «  "\n"; 

motr[j].set_initial(VM,  VA,  P,  Q,  X,  R,  1); 
eaf[j]  =  motr[j]  get_eaf(); 
tm[j]  =  -  motr[j].get_tm(), 
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cout  «  "Eaf  =  "  «  eaf[j]  «  "  Torque  =  "  «  tm[j]  «  "\n", 
} 

//Report  the  initial  conditions 

cout  «  "\n  Final  Setup  Report An", 
cout  «  "  Motor  An"; 

for  (j=0;  j<gen_count;  j++) 

{ 

motr[j].report(); 
cout  «  "\n"; 

} 
cout  «  "  Network  Nodes  An", 
for  (j=0;  j<net_node_count,  j++)  nptrfj]  report(); 
cout  «  "\n  Voltage  Source  Nodes:\n"; 
for  (j=0;  j<v_node_count;  j++)  vptr[j].report(); 
cout«"\nLines:\n"; 
for  (i=0;  i<nlines;  i++)  lptr[i].report(); 
cout «  "\n", 

//  Start  the  actual  time-step  simulation 

int  n,  m,  eptr=0,  flag=l,  ip=0, 

double  time=0, 

double  te[gen_count],  vbus[gen_count]; 

//  Output  first  data  point  (the  initial  conditions) 

for  (j=0;  j<gen_count,  j++)  //  Machine  output 

motr[j]  fileoutput(time), 

lineout  «  time  «  "  "; 

for  (i=0;  Knlines;  i++)  //  Output  to  line  file 

{ 
if  (line_out_key[i]  &  4)  lineout  «  lptr[i].i_a()  «  "  "; 

if  (line_out_key[i]  &  2)  lineout  «  lptr[i].i_b()  «  "  "; 

if  (line_out_key[i]  &  1)  lineout  «  lptr[i].i_c()  «  "  "; 

} 
lineout  «  "\n", 

nodeout  «  time  «  "  "; 

for  (j=0;  j<nnodes,  j++) 
switch(t[j]) 
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case  'g': 
if  (node 

«  "  "; 

out 

if  (node 

«  "  "; 

out 

if  (node 

«  "  "; 

out 

break; 

case  'n': 
if  (node 

«  "  "; 

out 

if  (node 

«  "  "; 

out 

if  (node 

«  "  "; 

out 

break, 

case  V: 
if  (node 

«  "  "; 

out 

if  (node 

«  "  "; 

out 

if  (node 

«  "  "; 

out 

break, 
default: 
break, 

} 
nodeout «  "\n"; 

keyjj]  &  4)  nodeout  «  motr[ptr[j]]  a_voltage() 
key[j]  &  2)  nodeout  «  motr[ptr[j]].b_voltage() 
key[j]  &  1)  nodeout  «  motr[ptr[j]]  cvoltageQ 

key[j]  &  4)  nodeout  «  nptr[ptr[j]]  a_voltage() 
key[j]  &  2)  nodeout  «  nptr[ptr[j]]  b_voltage() 
key[j]  &  1)  nodeout  «  nptr[ptr[j]].c_voltage() 

key[j]  &  4)  nodeout  «  vptr[ptr[j]]  a_voltage() 
key[j]  &  2)  nodeout  «  vptr[ptr[j]].b_voltage() 
keyO]  &  1)  nodeout  «  vptr[ptr[j]].c_voltage() 


double  err,  tol=le-8; 

for  (n=0;  n<no;  n++)  //  Outer  loop:  print  loop 

{ 
for  (m=0,  m<np;  m++) 

{ 
time  =  dt  *  (n  *  np  +  m);     //  Time  counter 

if  (time>=0.4)  { 

for  (j=0;j<gen_count,j++)  { 

tm[j]=  motr[j]  t_sc()*sqrt(0.4/time); 

} 


} 
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if  (time>=  event_time[eptr])     //  Event  monitor 
{ 
switch(event_type[eptr]) 

{ 
case  V: 

cout  «  "Opening  Line  "  «  event_line[eptr]  «  "  at  time  " 

« time  «  "\n"; 
lptr[event_line[eptr]].open_a(), 
lptr[event_line[eptr]]  open_b(); 
lptr[eventJine[eptr]].open_c(), 
break; 
case  'a': 
cout  «  "Closing  Phase  A  of  Line  "  «  eventline  [eptr] 

«  "  at  time  "  «  time  «  "\n", 
lptr[event_line[eptr]].close_a(); 
lptr[event_line[eptr]].setup(); 
flag=l; 
break; 
case  'b': 
cout  «  "Closing  Phase  B  of  Line  "  «  eventline  [eptr] 

«  "  at  time  "  «  time  «  "\n", 
lptr[event_line[eptr]]  close_b(), 
lptr[event_line[eptr]] .  setup(), 
flag=l; 
break; 
case  'c': 
cout  «  "Closing  Phase  C  of  Line  "  «  eventline  [eptr] 

«  "  at  time  "  «  time  «  "\n", 
lptr[event_line[eptr]]  close_c(), 
lptr[event_line[eptr]].setup(), 
flag=l; 
break; 
case  't': 
cout  «  "Closing  All  3  Phases  of  Line  "  «  eventline  [eptr] 

«  "  at  time  "  «  time  «  "\n", 
lptr[event_line[eptr]].close_a(), 
lptr[event_line[eptr]].close_b(); 
lptr[event_line[eptr]].close_c(); 
lptr[event_line[eptr]].setup(); 
flag=l; 
break, 

} 
eptr++, 

} 
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for  (i=0;  i<nlines;  i++) 
if  (ip  =  lptr[i].current_monitor()) 

{ 

lptr[i].re_setup(ip); 

flag=l; 
} 

if  (flag) 

{ 
for  (j=0;  j<net_node_count;  j++) 

nptr[j].setup(), 
flag  =  0; 

} 

if  (time>=0.4)  {  // Motor  windmilling 

omz  =  377*sqrt(0.4/time); 
} 

//  Runge-Kutta  step  1 

for  (j=0;  j<gen_count;  j++) 

motr[j].rkl(eaflj],  tm[j],  dt,  time), 
for  (j=0;  j<v_node_count,  j++) 

vptr[j].set_v(omz  *  time); 
err  =  1 ; 
do 

{err  =  0; 

for  (j=0;  j<net_node_count;  j++) 
err  +=  nptr[j]  estimatevoltage(dt), 

}  while  (err  >  tol); 

for  (i=0;  i<nlines,  i++) 
lptr[i].rkl(dt); 

//  Runge-Kutta  step  2 

time  =  dt*(n*np  +  m  + 0.5);     //Time 

for  (j=0;  j<gen_count,  j++) 
motr[j].rk2(eafTj],  tm[j],  dt,  time), 

for  (j=0;  j<v_node_count;  j++) 
vptr[j].set_v(omz  *  time), 

err  =  1 ; 
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do 

{err  =  0, 

for  (j=0;  j<net  nodecount;  j++) 
err  +=  nptr[j].estimate_voltage(dt); 

}  while  (err  >  tol); 
for  (i=0;  i<nlines,  i++) 

lptr[i].rk2(dt); 

//  Runge-Kutta  step  3 

for  (j=0;  j<gen_count,  j++) 

motr[j] .rk3(eaf[j],  tm[j],  dt,  time); 
err  =  1 , 
do 
{err  =  0; 
for  (j=0;  j<net_node_count;  j++) 

err  +=  nptr[j]  estimatevoltage(dt); 
}  while  (err  >  tol), 

for  (i=0;  i<nlines;  i++) 
lptr[i].rk3(dt); 

//  Runge-Kutta  step  4 

time  =  dt  *  (n  *  np  +  m  +  1 .0);     //  time 

for  (j=0;  j<gen_count,  j++) 

motr[j].rk4(eafU],  tm[j],  dt,  time), 
for  (j=0;  j<v_node_count;  j++) 

vptr[j]  set_v(omz  *  time); 
err  =  1 , 
do 

{err  =  0; 

for  (j=0,  j<net_node_count,  j++) 
err  +=  nptrQ]  estimatevoltage(dt); 

}  while  (err  >  tol), 

for  (i=0;  i<nlines;  i++) 
lptr[i].rk4(dt); 

//  Wrapup  the  Runge-Kutta  routine: 

for  (j=0;  j<gen_count;  j++) 

motr[j].rk(time); 
for  (i=0;  i<nlines,  i++) 
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lptr[i].rk(); 

//  Now  update  the  exciter  and  avr: 

for  (j=0;  j<gen_count;  j++) 
{ 

vbus[j]  =  motr[j].vsense(), 
te[j]  =  -motr|j].t_e(); 
} 
}  //  This  is  the  end  of  the  internal  loop 

for  (j=0;  j<gen_count;  j++)  //  Motor  output  here 

motr[j].file_output(time); 

lineout  «  time  «  "  "; 

for  (i=0;  i<nlines,  i++)  //  Output  to  line  file 

{ 

if  (line_out_key[i]  &  4)  lineout  «  lptr[i].i_a()  «  "  "; 

if  (line_out_key[i]  &  2)  lineout  «  lptr[i].i_b()  «  "  "; 

if  (line_out_key[i]  &  1)  lineout  «  lptr[i].i_c()  «  "  "; 

} 
lineout  «  "\n"; 

nodeout  «  time  «  "  "; 

for  (j=0,  j<nnodes,  j++) 
switch(t[j]) 

{ 

case  'g': 
if  (nodeoutkeyfj]  &  4)  nodeout  «  motr[ptr[j]].a_voltage() 


«  " 

«  " 

«  "  " 
break, 
case  'n': 


5 

f  (node_out_key[j]  &  2)  nodeout  «  motr[ptr[j]]  b_voltage() 
;f  (node_out_key[j]  &  1)  nodeout  «  motr[ptr[j]]  c_voltage() 


if  (node_out_key[j]  &  4)  nodeout  «  nptr[ptr[j]]  a_voltage() 
if  (node_out_key[j]  &  2)  nodeout  «  nptr[ptr[j]]  b_voltage() 
f  (node_out_key[j]  &  1)  nodeout  «  nptr[ptr[j]].c_voltage() 


«  "  "; 

break, 
case  V: 
if  (node_out_key[)]  &  4)  nodeout  «  vptr[ptr[j]].a_voltage() 
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«  "  "; 

if  (node_out_key[j]  &  2)  nodeout  «  vptr[ptr[j]].b_voltage() 
«  "  "; 

if  (node_out_key[j]  &  1)  nodeout  «  vptr[ptr[j]]  c_voltage() 

«  "  "■ 

break; 

default: 

break; 

} 


nodeout  «  "\n"; 


}  //  End  of  the  outer,  or  print  loop 

delete  [nlines]lptr; 
delete  [gen_count]motr, 
delete  [net_node_count]nptr, 
delete  [v_node_count]vptr; 
delete  busa, 
delete  busb, 

fn.close(), 
fl  close(), 


} 


void  concat  (char  *stringl,  char  *string2,  char  *string) 

{ 
for  (int  j=0;  j<strlen(stringl),  j++) 

string[j]  =  stringlO], 
for  (int  i=0;  i<strlen(string2),  i++) 
string[i+j]  =  string2[i]; 

} 
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Appendix  E-l.  Motor  Object 

The  following  program  models  the  permanent  magnet  motor    It  defines  the  class 
motor  that  is  executed  as  part  of  the  simulation  program    This  file  was  written  by  F  R 
Colberg  based  on  the  original  file  written  for  generators  by  Professor  James  L  Kirtley 

#include  <stream.h> 
#include  <math  h> 
#include  <libc.h> 
#include  "linen" 
#include  <stdlib.h> 


class  motor  :  public  bbus{ 
double  xd,  xq,  xdpp,  xqpp,  //  Reactances 

tdopp,  tqopp,  //  Time  Constants 

xz,  ta,  //  Armature  quantities 

h,  omz,  thz,  //  Miscellaneous  quantities 

eqpp,  edpp,  om,  delt,  //  State  Variables 

kqpp  1 ,  kdpp  1 ,  kdelt  1 ,  kom  1 ,  //  R-K  derivatives 
kqpp2,  kdpp2,  kdelt2,  kom2,  //  R-K  derivatives 
kqpp3,  kdpp3,  kdelt3,  kom3,  //  R-K  derivatives 
kqpp4,  kdpp4,  kdelt4,  kom4,        //  R-K  derivatives 


ea,  eb,  ec, 

ddl,  dd2,  dd3,  dd4, 

dql,  dq2,  dq3,  dq4, 

tmr,  eafr,tm, 

iaz,  ibz,  icz, 

tel,  te2,  te3,  te4,  te; 

int  node; 

line*  lptr; 

char  ofname[10]; 

filebuffl; 

public: 
void  genset(filebuf  fD) 


! 


istream  from  (&fD); 
from  »  xd; 
from  »  xq; 
from  »  xdpp; 
from  »  xqpp; 
from  »  xz; 
from  » tdopp, 
from  » tqopp; 


//  Output  Voltages 

//  intermediate  derivatives 

//  required  torque  and  eaf 
//  initial  currents 

//  electrical  torque 
//  node  number  of  internal  bus 
//  line  connecting  machine  to  net 
//  file  name  for  output 


//  Motor  reads  input  from  file 
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from  » ta, 
from  »  h; 
from  »  omz; 
from  »  ofhame; 
fl.open(ofname,  output); 
} 

void  setline  (line*  lineptr) 

{ 
lptr  =  lineptr; 

} 

void  setnode  (int  nodenumber) 

{ 
node  =  nodenumber; 

} 

void  file_output(double  time) 

{ 

ostream  to  (&fl); 

to  «  time  «  "  "  «  eqpp  «  "  "  «  edpp  «  " 
«  "  "  «  delt  «  "  "  «  om  «  "\n", 

} 

void  report() 

{ 
cout «  "  MotorVn" 

«  "  xd  =  "  «  xd 

«  "  xq  =  "  «  xq 

«  "  xdpp  =  "  «  xdpp 

«  "  xqpp  =  "  «  xqpp  «  "\n" 

«  "  tdopp  =  "  «  tdopp 

«  "  tqopp  =  "  «  tqopp 

«  "  h  =  "  «  h 

«  "\n" 

«  "  eqpp  =  "  «  eqpp 

«  "  edpp  =  "  «  edpp 

«  "  om  =  "  «  om 

«  "  delt  =  "  «  delt  «  "V 


} 


«  "  node  =  "  «  node 

«  "  conn  to  bus  "  «  lptr->report_node_b()  «  "\n", 


void  set  initial  (double  vl,  double  thl,  double  pi,  double  ql, 
double  xl,  double  rl,  int  longout  =  0) 
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//  set  initial  conditions  from 
//  external  bus.  xl  and  rl  are 
//  branch  inductance  (including 
//  subtransient  reactance ) 

{ 

double  psi  =  atan2  (ql,  pi);        //  power  factor  angle  at  ext  bus 
double  il  =  sqrt(pl*pl+ql*ql)/vl;  //  load  current 
double  xt  =  xq  +  xl  -  xdpp,         //  total  inductance  to  ext  bus 
double  vr  =  vl  -  rl  *  il  *  cos  (psi)  +  xt  *  il  *  sin  (psi), 
double  vi  =  xt  *  il  *  cos  (psi)  +  rl  *  il  *  sin  (psi); 
double  vs  =  sqrt  (vr*vr  +  vi*vi); 
double  ths  =  thl  -  atan2  (vi,  vr), 
double  th,  ca,  cb,  cc,  sa,  sb,  sc; 
doubles  =  2.0943951;        // 2  pi/3 

thz  =  thl-  1.5707963; 
delt  =  ths  -  thl; 

double  id  =  il  *  sin  (delt  +  psi);         //  direct  and  quadrature  axis  currents 
double  iq  =  il  *  cos  (delt  +  psi); 
eafr  =  vs  +  id  *  (xd  -  xq); 
om  =  omz; 

eqpp  =  eafr  +  (xd  -  xdpp)  *  id; 
edpp  =  -(xq  -  xqpp)  *  iq; 

tmr  =  (eqpp*iq  +  edpp*id  +  (xdpp  -  xqpp)  *  id  *  iq), 

iaz  =  il  *  cos  (thl  +  psi); 

ibz  =  il  *  cos  (thl  +  psi  -  2.0943951), 

icz  =  il  *  cos  (thl  +  psi  +  2.0943951); 

th  =  thz  +  delt; 

ca  =  cos  (th);  cb  =  cos  (th  -  s);  cc  =  cos  (th  +  s); 

sa  =  sin  (th);  sb  =  sin  (th  -  s);  sc  =  sin  (th  +  s); 

ea  =  -  (om/omz)*(eqpp*sa-edpp*ca), 
eb  =  -  (om/omz)*(eqpp*sb-edpp*cb), 
ec  =  -  (om/omz)*(eqpp*sc-edpp*cc), 

set_va(ea);  set_vb(eb);  setvc(ec); 

if  (longout) 

{ 

cout  «  form("mot::set_initial()  Node  =  %d\n",  node); 
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cout  «  form("vl  =  %g  thl  =  %g  pi  =  %g  ql  =  %g\n",  vl,  thl,  pi,  ql), 
cout  «  form("xl  =  %g  rl  =  %g    psi  =  %g  il  =  %g\n",  xl,  rl,  psi,  il), 
cout  «  form("xt  =  %g  vr  =  %g  vi  =  %g  vs  =  %g\n",  xt,  vr,  vi,  vs), 
cout  «  form("ths  =  %g  thz  =  %g  delt  =  %g  id  =  %g\n",ths,thz,delt,  id), 
cout  «  form("iq  =  %g  eafr  =  %g  tmr  =  %g  om  =  %g\n",  iq,eafr,  tmr,  om), 
cout  «  form("eqpp  =  %g  edpp  =  %g\n\n",  eqpp,  edpp), 
} 


} 


double  get_tm()  {return  (tmr);} 
double  get_eaf()  {return  (eafr),} 
double  get_ia()  {return  (iaz),} 
double  get_ib()  {return  (ibz),} 
double  get_ic()  {return  (icz);} 
double  get_eqpp()  {return  (eqpp);} 
double  get_edpp()  {return  (edpp),} 
double  t_e()  {return  (te);} 
double  t_el()  {return  (tel);} 
double  t_e2()  {return  (te2);} 
double  t_e3()  {return  (te3),} 
double  t_e4()  {return  (te4);} 

double  t_sc()  {  //Calculates  input  torque  to  the 

tm  =  get_tm();  //motor  with  the  shaft  free  wheeling 

return(tm); 
} 

double  eaf_sc()  { 
double  eaf  =  get_eaf(); 
return  (eaf); 
} 

void  disp_params()  {cout  «  "  xd  =  "  «  xd 

«  "  xq  =  "  «  xq 

«  "\n" 

«  "  xdpp  =  "  «  xdpp 

«  "  xqpp  =  "  «  xqpp  «  "\n" 

«  "  tdopp  =  "  «  tdopp 

«  "  tqopp  =  "  «  tqopp  «  V 

«  "  h  =  "  «  h 

«  "  omz  =  "  «  omz  «  "  thz  ="  «  thz  «"\n",} 
void  disp_state()   {cout  «  "eqpp  =  "  «  eqpp 

«  "\tedpp  =  "  «  edpp  «  "\n" 

«  "  om  =  "  «  om  «  "\tdelta  =  "  «  delt  «  "\n";  } 
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double  e_a()  (return(ea);} 
double  e_b()  (return(eb),} 
double  e_c()  (return(ec),} 
double  omega()  { return(om), } 
double  delta()  {return(delt);} 

void  rkl  (double  eaf,  double  tm, 

double  dt,  double  t,  int  longout  =  0) 

{ 

double  ia,  ib,  ic, 
double  th,  id,  iq,  deqpp,  dedpp,  ddelt,  dom, 

ca,  cb,  cc,  sa,  sb,  sc, 
double  s  =  2.094395 1 ,        //  2  pi  /  3 

//  First,  Park's  Transform  Currents 

ia  =  -lptr->i_ap(); 
ib  =  -lptr->i_bp(); 
ic  =  -lptr->i_cp(); 

if  (longout) 

{ 

cout  «form("rk  l:\n"); 
cout «  form("node  %d\n",node); 
cout  «  formfia  =  %g  ib  =  %g  ic  =  %g\n",  ia,  ib,  ic), 
cout  «  form("eaf  =  %g  tm  =  %g  dt  =  %g  t  =  %g\n",  eaf,tm,dt,t); 

} 

th  =  thz  +  omz  *  t  +  delt; 

ca  =  cos  (th);  cb  =  cos  (th  -  s);  cc  =  cos  (th  +  s), 

sa  =  sin  (th),  sb  =  sin  (th  -  s);  sc  =  sin  (th  +  s); 

id  =  6666666667  *  ( ia  *  ca  +  ib  *  cb  +  ic  *  cc  ), 
iq  =  6666666667  *  ( -ia  *  sa  -  ib  *  sb  -  ic  *  sc); 

if  (longout) 

{ 

cout  «  form("th  =  %g  ca  =  %g  sa  =  %g\n",  th,  ca,  sa), 

cout  «  formfid  =  %g  iq  =  %g\n",  id,  iq); 
} 

//  Now  do  the  time  step 
deqpp  =  -  eqpp/tdopp  +  eaf/tdopp  +  (xd  -  xdpp)  *  id  /  tdopp, 
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dedpp  =  -  edpp/tqopp  -  (xq  -  xqpp)  *  iq  /  tqopp, 
ddelt  =  om  -  omz; 
if  (t>=0.4)  { 
tel  =0.0; 

} 
else 

{tel  =  eqpp*iq  +  edpp*id  +(xdpp  -  xqpp)  *  id  *  iq;} 
dom  =  (omz/(2.0*h))  *  (tel  +  tm), 

if(longout) 

{ 

cout  «  form("deqpp  =  %g  dedpp  =  %g\nM,  deqpp,  dedpp), 
cout  «  form( "ddelt  =  %g  dom  =  %g\n",  ddelt,  dom); 

} 

//  Now  the  R-K  coefficients  are: 

kqppl  =dt  *  deqpp, 
kdppl  =  dt  *  dedpp; 
kdeltl  =  dt  *  ddelt; 
koml  =  dt  *  dom; 

//  save  for  final  step 

dq 1  =  deqpp; 
ddl  =  dedpp, 

if  (long  out) 

{ 

cout  «  form("kqpp  =  %g  kdpp  =  %g\n",  kqppl,  kdppl), 
cout  «  form("kdelt  =  %g  kom  =  %g\n",  kdeltl,  koml), 

} 

ea  =  -  (om/omz)*((eqpp+.5*kqppl)*sa-(edpp+  5*kdppl)*ca) 

+  (1  0/omz)*(deqpp*ca+dedpp*sa); 
eb  =  -  (om/omz)*((eqpp+.5*kqppl)*sb-(edpp+.5*kdppl)*cb) 

+  (1.0/omz)*(deqpp*cb+dedpp*sb), 
ec  =  -  (om/omz)*((eqpp+.5*kqppl)*sc-(edpp+.5*kdppl)*cc) 

+  (1  0/omz)*(deqpp*cc+dedpp*sc); 

setva(ea), 
setvb(eb); 
setvc(ec), 
if  (t<0.3)  { 
set_varc(0.5*((ea-eb)));} 
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if(longout) 

{ 
cout  «  form("ea  =  %g  eb  =  %g  ec  =  %g\n\n",  ea,  eb,  ec); 

} 
} 

void  rk2  (double  eaf,  double  tm, 

double  dt,  double  t,  int  longout  =  0) 

{ 
double  th,  id,  iq,  deqpp,  dedpp,  ddelt,  dom, 

ca,  cb,  cc,  sa,  sb,  sc; 
double  s  =  2.094395 1 ;        //  2  pi  /  3 
double  eqppt,  edppt,  deltt,  omt,    //  temporaries  for  states 
double  ia,  ib,  ic, 

//  assign  state  temporaries: 

eqppt  =  eqpp  +  .  5  *  kqpp  1 ; 
edppt  =  edpp  +  5  *  kdppl; 
omt  =  om  +  .5  *  koml; 
deltt  =  delt  +  .5  *  kdeltl; 

//  First,  Park's  Transform  Currents 

ia  =  -lptr->i_ap(); 
ib  =  -lptr->i_bp(); 
ic  =  -lptr->i_cp(), 

if  (longout) 

{ 

cout «  form("rk2:\n"); 

cout  «  form("node  %d\n",  node), 

cout  «  formfia  =  %g  ib  =  %g  ic  =  %g\n",  ia,  ib,  ic), 

cout  «  formfeaf  =  %g  tm  =  %g  dt  =  %g  t  =  %g\n",  eaf,tm,dt,t), 
} 

th  =  thz  +  omz  *  t  +  delt; 

ca  =  cos  (th);  cb  =  cos  (th  -  s),  cc  =  cos  (th  +  s); 

sa  =  sin  (th),  sb  =  sin  (th  -  s),  sc  =  sin  (th  +  s); 

id  =  6666666667  *  ( ia  *  ca  +  ib  *  cb  +  ic  *  cc  ), 
iq  =  .6666666667  *  ( -ia  *  sa  -  ib  *  sb  -  ic  *  sc), 
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if(longout) 

{ 

cout  «  form("th  =  %g  ca  =  %g  sa  =  %g\n",  th,  ca,  sa); 

cout  «  form("id  =  %g  iq  =  %g\n",  id,iq); 
} 

//  Now  do  the  time  step 

deqpp  =  -  eqppt/tdopp  +  eaf/tdopp  +  (xd  -  xdpp)  *  id  /  tdopp, 

dedpp  =  -  edppt/tqopp  -  (xq  -  xqpp)  *  iq  /  tqopp, 

ddelt  =  omt  -  omz, 

if  (t>=0  4)  { 

te2=0.0;} 

else 

{te2  =  eqppt*iq  +  edppt*id  +(xdpp  -  xqpp)  *  id  *  iq;} 

dom  =  (omz/(2.0*h))  *  (te2  +  tm), 

if  (longout) 

{ 

cout  «  form("deqpp  =  %g  dedpp  =  %g\n",  deqpp,  dedpp); 

cout  «  form("ddelt  =  %g  dom  =  %g\n",  ddelt,dom), 
} 

//  Now  the  R-K  coefficients  are: 

kqpp2  =  dt  *  deqpp, 
kdpp2  =  dt  *  dedpp, 
kdelt2  =  dt  *  ddelt, 
kom2  =  dt  *  dom; 

//  save  for  final  step 

dq2  =  deqpp; 
dd2  =  dedpp; 

if  (longout) 

{ 

cout  «  form("kqpp  =  %g  kdpp  =  %g\n",  kqppl,  kdppl), 
cout  «  form("kdelt  =  %g  kom  =  %g\n",  kdeltl,  koml); 

} 

ea  =  -  (om/omz)*((eqpp+.5*kqpp2)*sa-(edpp+.5*kdpp2)*ca) 

+  ( 1  0/omz)*(deqpp*ca+dedpp*sa); 
eb  =  -  (om/omz)*((eqpp+.5*kqpp2)*sb-(edpp+.5*kdpp2)*cb) 

+  (1.0/omz)*(deqpp*cb+dedpp*sb); 
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ec  =  -  (om/omz)*((eqpp+.5*kqpp2)*sc-(edpp+.5*kdpp2)*cc) 
+  (1  0/omz)*(deqpp*cc+dedpp*sc); 

set_va(ea), 
set_vb(eb), 
set_vc(ec), 

if  (t<0,3)  { 
set_varc(0. 5  *((ea-eb)), } 

if  (longout) 

{ 
cout  «  form("ea  =  %g  eb  =  %g  ec  =  %g\n\n",  ea,  eb,  ec); 

} 
} 

void  rk3  (double  eaf,  double  tm, 

double  dt,  double  t,  int  longout  =  0) 

{ 
double  th,  id,  iq,  deqpp,  dedpp,  ddelt,  dom, 

ca,  cb,  cc,  sa,  sb,  sc, 
doubles  =  2.0943951;        // 2  pi  /  3 
double  eqppt,  edppt,  deltt,  omt;  //  temporaries  for  states 
double  ia,  ib,  ic, 

//  assign  state  temporaries: 

eqppt  =  eqpp  +  .5  *  kqpp2, 
edppt  =  edpp  +  .5  *  kdpp2; 
omt  =  om  +  .5  *  kom2; 
deltt  =  delt  +  5  *  kdelt2; 

//  First,  Park's  Transform  Currents 

ia  =  -lptr->i_ap(), 
ib  =  -lptr->i_bp(), 
ic  =  -lptr->i_cp(), 

if  (longout) 

{ 

cout«form("rk3:\n"), 
cout  «  form("node  %d\n",  node), 
cout  «  form("ia  =  %g  ib  =  %g  ic  =  %g\n",  ia,  ib,  ic), 
cout  «  formfeaf  =  %g  tm  =  %g  dt  =  %g  t  =  %g\n",  eaf,tm,dt,t), 
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} 

th  =  thz  +  omz  *  t  +  delt, 

ca  =  cos  (th);  cb  =  cos  (th  -  s);  cc  =  cos  (th  +  s); 

sa  =  sin  (th);  sb  =  sin  (th  -  s);  sc  =  sin  (th  +  s), 

id  =  .6666666667  *  ( ia  *  ca  +  ib  *  cb  +  ic  *  cc  ); 
iq  =  .6666666667  *  ( -ia  *  sa  -  ib  *  sb  -  ic  *  sc), 

if  (longout) 

I 

cout  «  form("th  =  %g  ca  =  %g  sa  =  %g\n",  th,  ca,  sa); 

cout  «  form("id  =  %g  iq  =  %g\n",  id,iq); 
} 

//  Now  do  the  time  step 

deqpp  =  -  eqppt/tdopp  +  eaf/tdopp  +  (xd  -  xdpp)  *  id  /  tdopp, 
dedpp  =  -  edppt/tqopp  -  (xq  -  xqpp)  *  iq  /  tqopp, 
ddelt  =  omt  -  omz; 

if  (t>=0  4)  { 

te2=0.0;} 
else 

{te3  =  eqppt*iq  +  edppt*id  +(xdpp  -  xqpp)  *  id  *  iq;} 
dom  =  (omz/(2.0*h))  *  (te3  +  tm); 

//  save  for  final  step 

dq3  =  deqpp, 
dd3  =  dedpp; 

if  (longout) 

{ 

cout  «  form("deqpp  =  %g  dedpp  =  %g\n",  deqpp,  dedpp); 

cout  «  form("ddelt  =  %g  dom  =  %g\n",  ddelt,dom), 
} 

//  Now  the  R-K  coefficients  are: 

kqpp3  =  dt  *  deqpp, 
kdpp3  =  dt  *  dedpp; 
kdelt3  =  dt  *  ddelt; 
kom3  =  dt  *  dom, 

if  (longout) 
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{ 

cout  «  form("kqpp  =  %g  kdpp  =  %g\n",  kqppl,  kdppl); 
cout  «  formfkdelt  =  %g  kom  =  %g\n",  kdeltl,  koml), 
} 

//  And  now  for  the  voltages: 


ea  =  -  (om/omz)*((eqpp+kqpp3)*sa-(edpp+kdpp3)*ca) 

+  (1.0/omz)*(deqpp*ca+dedpp*sa), 
eb  =  -  (om/omz)*((eqpp+kqpp3)*sb-(edpp+kdpp3)*cb) 

+  (1  0/omz)*(deqpp*cb+dedpp*sb); 
ec  =  -  (om/omz)*((eqpp+kqpp3)*sc-(edpp+kdpp3)*cc) 

+  (1.0/omz)*(deqpp*cc+dedpp*sc); 

set_va(ea), 
set_vb(eb), 
setvc(ec), 

if  (t<0.3)  { 
set_varc(0.5*((ea-eb)));} 

if  (longout) 

{ 
cout  «  form("ea  =  %g  eb  =  %g  ec  =  %g\n\n",  ea,  eb,  ec); 

} 


! 


void  rk4  (double  eaf,  double  tm, 

double  dt,  double  t,  int  longout  =  0) 

{ 
double  th,  id,  iq,  deqpp,  dedpp,  ddelt,  dom, 

ca,  cb,  cc,  sa,  sb,  sc; 
double  s  =  2.094395 1 ;        //  2  pi  /  3 
double  eqppt,  edppt,  deltt,  omt,  //  temporaries  for  states 
double  ia,  ib,  ic, 

//  assign  state  temporaries: 

eqppt  =  eqpp  +  kqpp3; 
edppt  =  edpp  +  kdpp3; 
omt  =  om  +  kom3, 
deltt  =  delt  +  kdelt3: 
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//  First,  Park's  Transform  Currents 

ia  =  -lptr->i_ap(); 
ib  =  -lptr->i_bp(); 
ic  =  -lptr->i  cp(); 

if  (longout) 

{ 

cout «  form("rk4:\n"), 

cout  «  form("node  %d\n",  node), 

cout  «  formfia  =  %g  ib  =  %g  ic  =  %g\n",  ia,  ib,  ic), 

cout  «  formfeaf  =  %g  tm  =  %g  dt  =  %g  t  =  %g\n",  eaf,tm,dt,t), 
} 

th  =  thz  +  omz  *  t  +  delt, 

ca  =  cos  (th);  cb  =  cos  (th  -  s);  cc  =  cos  (th  +  s); 

sa  =  sin  (th);  sb  =  sin  (th  -  s);  sc  =  sin  (th  +  s); 

id  =  .6666666667  *  ( ia  *  ca  +  ib  *  cb  +  ic  *  cc  ); 
iq  =  .6666666667  *  ( -ia  *  sa  -  ib  *  sb  -  ic  *  sc), 

if  (longout) 
{  cout  «  form("th  =  %g  ca  =  %g  sa  =  %g\n",  th,  ca,  sa); 

cout  «  formfid  =  %g  iq  =  %g\n",  id,iq), 
} 

//  Now  do  the  time  step 

deqpp  =  -  eqppt/tdopp  +  eaf/tdopp  +  (xd  -  xdpp)  *  id  /  tdopp, 
dedpp  =  -  edppt/tqopp  -  (xq  -  xqpp)  *  iq  /  tqopp, 
ddelt  =  omt  -  omz; 

if  (t>=0  4)  { 

te2=0.0;} 

else 

{te4  =  eqppt*iq  +  edppt*id  +(xdpp  -  xqpp)  *  id  *  iq,} 

dom  =  (omz/(2.0*h))  *  (te4  +  tm), 

if  (longout) 

{ 

cout  «  form("deqpp  =  %g  dedpp  =  %g\n",  deqpp,  dedpp), 
cout  «  form("ddelt  =  %g  dom  =  %g\n",  ddelt,dom), 

} 

//  The  R-K  coefficients  are: 
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kqpp4  =  dt  *  deqpp; 
kdpp4  =  dt  *  dedpp; 
kdelt4  =  dt  *  ddelt; 
kom4  =  dt  *  dom, 

if(longout) 

{ 

cout  «  form("kqpp  =  %g  kdpp  =  %g\n",  kqppl,  kdppl), 
cout  «  form("kdelt  =  %g  kom  =  %g\n",  kdeltl,  koml), 

} 

//  save  for  final  step 

dq4  =  deqpp; 
dd4  =  dedpp, 

ea  =  -  (om/omz)*(eqpp*sa-edpp*ca)+(l  0/omz)* (deqpp *ca+dedpp*sa); 
eb  =  -  (om/omz)*(eqpp*sb-edpp*cb)+(l  0/omz)*(deqpp*cb+dedpp*sb), 
ec  =  -  (om/omz)*(eqpp*sc-edpp*cc)+(l  0/omz)*(deqpp*cc+dedpp*sc); 

setva(ea), 
setvb(eb); 
setvc(ec), 

if  (t<0.3)  { 

set_varc(0.5*((ea-eb))),} 

if  (longout) 

{ 
cout  «  form("ea  =  %g  eb  =  %g  ec  =  %g\n\nM,  ea,  eb,  ec); 

} 
} 
//  End  the  Runge-Kutta  routine 

void  rk(double  t,  int  longout  =  0) 

{ 

double  th,  ca,  cb,  cc,  sa,  sb,  sc, 
double  ia,  ib,  ic,  id,  iq; 
double  s  =  2.094395 1 ,        //  2  pi  /  3 

ia  =  -lptr->i_ap(); 
ib  =  -lptr->i_bp(); 
ic  =  -lptr->i_cp(); 
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th  =  thz  +  omz  *  t  +  delt, 

ca  =  cos  (th);  cb  =  cos  (th  -  s),  cc  =  cos  (th  +  s); 

sa  =  sin  (th),  sb  =  sin  (th  -  s);  sc  =  sin  (th  +  s); 

id  =  6666666667  *  ( ia  *  ca  +  ib  *  cb  +  ic  *  cc  ); 
iq  =  .6666666667  *  (  -ia  *  sa  -  ib  *  sb  -  ic  *  sc), 

eqpp  +=  (kqppl  +  2*kqpp2  +  2*kqpp3  +  kqpp4)/6.0, 
edpp  +=  (kdppl  +  2*kdpp2  +  2*kdpp3  +  kdpp4)/6.0; 
delt  +=  (kdeltl  +  2*kdelt2  +  2*kdelt3  +  kdelt4)/6 .0; 
om  +=  (koml  +  2*kom2  +  2*kom3  +  kom4)/6.0; 

double  deqpp  =  (dql  +  2*dq2  +  2*dq3  +  dq4)/6, 
double  dedpp  =  (ddl  +  2*dd2  +  2*dd3  +  dd4)/6; 

ea  =  -  (orn/omz)*(eqpp*sa-edpp*ca)+(1.0/omz)*(deqpp*ca+dedpp*sa), 
eb  =  -  (om/omz)*(eqpp*sb-edpp*cb)+(l  0/omz)*(deqpp*cb+dedpp*sb); 
ec  =  -  (om/ornz)*(eqpp*sc-edpp*cc)+(1.0/omz)*(deqpp*cc+dedpp*sc), 

set_va(ea); 
set_vb(eb), 
setvc(ec), 

if  (t<0  3)  { 

set_varc(0.5*((ea-eb)));} 

te  =  (tel  +  2.0  *  te2  +  2.0  *  te3  +  te4)/6.0; 

if  (longout) 

{ 

cout  «  form("rk:\n"); 

cout  «  form("eqpp  =  %g  edpp  =  %g\n",  eqpp,  edpp), 

cout  «  form("delt  =  %g  om  =  %g  \n\n",  delt,  om), 
} 

filebuffr; 

char  int_volt[15]  -'int.dat", 
double  isc,  e_sc,  psc; 
fr.open(int_volt,  append), 
ostream  voltout(&fr); 
voltout«t«"  "; 

if((t>0.3)&&((0  5*(ea-eb)>0  05)||(0  5*(ea-eb)<-0  05)))  { 
i_sc  =  -0.5*(ea-eb)/xdpp*cos(om*t  +  acos(0  07854/(varc))) 

+  0.05/xdpp(l  57079-om*t);}       //Calculation  of  short  circuit  current 
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else  { 
i_sc  =  0.0;} 
if(i_sc>0.0) 

{e_sc=0.05;} 
else 
if  ( i_sc<  0.0) 

{e_sc=-0.05;} 
else 
{e_sc=  0.0} 


p_sc=e_sc*i_sc, 

voltout«e_sc«"  "«i  sc«"  ", 
voltout«p_sc«"  "; 
voltout«"\n" 

} 

}; 

//End  of  class  motor 
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Appendix  E-2.  Network  Program 

This  program  calculates  currents  and  voltages  of  network  buses  during  the 
simulation    It  is  executed  with  the  network  simulation  program,  Appendix  E. 

This  program  was  written  by  Professor  James  L  Kirtley 

^include  <stream.h> 
#include  <math  h> 

class  line; 

class  bbus  //  Base  class  for  polyphase  buses 

{ 
double  va,  vb,  vc,        //  All  buses  have  voltages  in  common 

public: 

double  a_voltage()      {  return(va);} 
double  b_voltage()      {  return(vb);} 
double  c_voltage()      {  return(vc);} 
void  set_va(double  v)  {va  =  v;} 
void  set_vb(double  v)  {vb  =  v;} 
void  set_vc(double  v)  {vc  =  v,} 
void  set_varc(double  v)  {return(va-vb);} 

}; 

class  bus  :  public  bbus  {        //  network  bus  for  time-step  simulations 
line**  lptr,  //  bus  is  connected  to  a  bunch  of  lines 

int  nodeno;  //  we  gotta  know  which  node  we  are 

int  nlines;  //  and  this  is  the  number  of  lines 

double  ia,  ib,  ic;  //  unbalanced  bus  currents 

double  gaai,  gabi,  gaci;       //  inverse  admittances 
double  gbai,  gbbi,  gbci, 
double  gcai,  gcbi,  gcci, 
double  gaat,  gabt,  gact; 
double  gbat,  gbbt,  gbct; 
double  gcat,  gcbt,  gcct, 

public: 

void  setbus(int  nodenumber,  int  linecount,  line*  linebuf[]); 
void  setup  ();  //  this  step  sets  up  node  admittances 

double  estimate_voltage(double  dt);    //  to  get  voltage  of  an  isolated  node 
void  reportQ, 
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};  //  end  of  declaration  of  class  bus 

class  vbus  :  public  bbus  {       //  voltage-source  bus 
double  vamp,  vphase, 

public: 

void  setvbus(double  v,  double  p) 

{ 

vamp  =  v; 

vphase  =  p; 

} 
void  report(); 
void  set_v(double  omt) 

{ 

double  s=2. 0943 951, 
set_va(vamp  *  cos  (omt  +  vphase)), 
set_vb(vamp  *  cos  (omt  +  vphase  -  s)); 
set_vc(vamp  *  cos  (omt  +  vphase  +  s)); 

} 
};  // end  of  declaration  of  vbus 

class  line  { 
double  xs,  xm,  iao,  ibo,  ico,  ia,  ib,  ic, 

omz,  kla,  k2a,  k3a,  k4a,  klb,  k2b,  k3b,  k4b, 

klc,  k2c,  k3c,  k4c,  iap,  ibp,  icp; 
int  nodea,  nodeb,  sa,  sb,  sc, 
bbus*  busa; 
bbus*  busb, 

public: 

double      gaa,  gab,  gac,  gba,  gbb,  gbc,  gca,  gcb,  gcc,  r; 

void  setline  (int  na,  int  nb,  double  x,  double  xz,  double  ra, 
double  omza,  int  saa,  int  sba,  int  sea, 
double  iaa,  double  iba,  double  ica); 

void  set_bus_pointers(bbus*  busa,  bbus*  bus_b), 

void  initcurrents  (double  iaa,  double  iba,  double  ica); 

void  setup  (),  //  Conductance  Parameters 

void  open_a(), 

void  open_b(), 

void  open_c(); 

void  close_a(); 

void  close_b(); 

void  close_c(), 

void  report(), 

void  report(char  *label), 
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int  currentmonitor  ();  //  Checks  for  zero  crossings 

void  resetup  (int  ip),    //  re-build  line  model 

void  rkl  (double  dt); 

void  rk2  (double  dt); 

void  rk3  (double  dt); 

void  rk4  (double  dt); 

void  rk  ();  //  Finishes  the  Runge-Kutta; 

double  i_a  ()  {  return  (ia);} 
double  i_b  ()  {  return  (ib);} 
double  i_c  ()  {  return  (ic),} 

double  ias(int  nn);  //  current  with  proper  sign  convention 

double  ibs(int  nn); 
double  ics(int  nn); 

double  i_ap()  (return(iap),}  //  Partial  deltas  for  runge-kutta  step 
double  i_bp()  (return(ibp);} 
double  i_cp()  {return(icp);} 

double  g_aa()  {return  (gaa), 
double  g_ab()  {return  (gab): 
double  g_ac()  {return  (gac); 
double  g_ba()  {return  (gba); 
double  g_bb()  {return  (gbb) 
double  g_bc()  {return  (gbc); 
double  g_ca()  {return  (gca); 
double  g_cb()  {return  (gcb) 
double  g_cc()  {return  (gcc): 

int  report_node_a(); 
int  report_node_b(), 

double  vaotherend  (int  thisnode); 
double  vbotherend  (int  thisnode); 
double  vcotherend  (int  thisnode); 

int  reportothernode  (int  thisnode); 

int  abs(int  x)  {  if(x<0)  return(-x);  else  return(x);} 

};  // end  of  definition  of  class  line 

void  bus::setbus(int  nodenumber,  int  linecount,  line*  linebuf[]) 
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} 


nodeno  =  node  number; 
nlines  =  linecount; 
lptr  =  new  line*[line_count], 
for  (int  j=0;  j<nlines,  j++) 

lptr[j]  =  linebuflj]; 
ia=ib=ic=0: 


void  bus:  setup  () 

{ 

gaat=0;  gabt=0,  gact=0; 
gbat=0,  gbbt=0,  gbct=0, 
gcat=0;  gcbt=0;  gcct=0; 


//  this  step  sets  up  node  admittances 


for  (int  i= 

{ 

gaat  +: 
gabt  +-- 
gact  +: 
gbat  +: 
gbbt  + 
gbct  +: 
gcat  +: 
gcbt  +: 
gcct  +: 

} 


0;  i<nlines,  i++) 


lptrfi]- 
lptr[i]- 
lptr[i]- 
lptr[i]- 
lptr[i] 
lptr[i] 
lptr[i]- 
lptr[i] 
lptr[i]- 


>gaa, 
>gab; 
>gac, 
>gba, 
■>gbb, 
>gbc; 
>gca; 
•>gcb; 
>gcc; 


//  invert  that:  since  it  is  a  3x3,  we  do  directly 

double  det  =  gaat* gbbt* gcct  +  gabt* gbct* gcat  +  gact* gbat* gcbt 
-  gaat*gcbt*gbct  -  gbat* gabt* gcct  -  gcat* gbbt* gact, 


gaai 
gabi 
gaci 
gbai 
gbbi 
gbci 
gcai 
gcbi 
gcci 


(gbbt*gcct 
(gcat*  gbct 
(gbat*  gcbt 
(gcbt* gact 
(gaat*  gcct 
(gcat*  gabt 
(gabt*  gbct 
(gbat*gact 
(gaat*  gbbt 


gcbt*gbct)/det; 
gbat*gcct)/det, 
gcat*gbbt)/det, 
gabt*gcct)/det; 
gcat*gact)/det, 
gaat*  gcbt  )/det, 
gbbt*gact)/det, 
gaat*  gbct  )/det; 
gbat*gabt)/det, 


cout  «  "bus::setup()  Here  is  Total  Admittance  G_t\n"; 
cout  «  gaat  «  "  "  «  gabt  «  "  "  «  gact  «  V; 
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cout  «  gabt  «  "  "  «  gbbt  «  "  "  «  gbct  «  "\n"; 
cout  «  gcat  «  "  "  «  gcbt  «  "  "  «  gcct  «  "\n", 

cout  «  "And  inverse  GOAn"; 
cout  «  gaai  «  "  "  «  gabi  «  "  "  «  gaci  «  "\n"; 
cout  «  gbai  «  "  "  «  gbbi  «  "  "  «  gbci  «  "\n", 
cout  «  gcai  «  "  "  «  gcbi  «  "  "  «  gcci  «  "\n"; 
} 

double  bus::estimate_voltage(double  dt)  //  voltage  of  an  isolated  node 

{ 

double  gva, 
double  gvb, 
double  gvc; 
int  i; 

double  ia,  ib,  ic; 

double  iau,  ibu,  icu,      //  estimated  unbalance  currents 
double  ova,  ovb,  ovc; 

double  va,  vb,  vc;  //  line  active  voltage 

double  vca,  vcb,  vcc,         //  correction  voltages 
double  err; 

iau=ibu=icu=gva=gvb=gvc=0 ; 

ova  =  a_voltage(), 
ovb  =  b_voltage(); 
ovc  =  c_voltage(); 
//        cout  «  "bus::estimate_voltage("«  nodeno«  ")\n"; 

for  (i=0;  i<nlines;  i++) 

{ 

ia  =  lptr[i]->ias(nodeno); 
ib  =  lptr[i]->ibs(nodeno); 
ic  =  lptr[i]->ics(nodeno); 

va  =  lptr[i]->va_other_end(nodeno)  +  lptr[i]->r  *  ia; 
vb  =  lptr[i]->vb_other_end(nodeno)  +  lptr[i]->r  *  ib; 
vc  =  lptr[i]->vc_other_end(nodeno)  +  lptr[i]->r  *  ic; 

//  cout  «  "Node  "  «  lptr[i]->report_other_node(nodeno) 

//  «  -  v  =  "  «  lptr[i]->va_other_end(nodeno)  « 

//  «  lptr[i]->vb_other_end(nodeno)  «  "  " 

//  «  lptr[i]->vc_other_end(nodeno)  «  "\n* 


M  ii 


gva  +=  lptr[i]->gaa  *  va 
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+  lptr[i]->gab  *  vb 
+  lptr[i]->gac  *  vc; 

gvb  +=  lptr[i]->gba  *  va 
+  lptr[i]->gbb  *  vb 
+  lptr[i]->gbc  *  vc; 

gvc  +=  lptr[i]->gca  *  va 
+  lptr[i]->gcb  *  vb 
+  lptr[i]->gcc  *  vc; 

iau  +=  ia, 
ibu  +=  ib; 
icu  +=  ic; 

} 

double  vcac,  vcbc,  vccc; 
double  vaf,  vbf,  vcf; 

vca  =  gaai*gva  +  gabi*gvb  +  gaci*gvc, 
vcb  =  gbai*gva  +  gbbi*gvb  +  gbci*gvc, 
vcc  =  gcai*gva  +  gcbi*gvb  +  gcci*gvc, 

vcac  =  -  (5/dt)*(gaai*iau  +  gabi*ibu  +  gaci*icu), 
vcbc  =  -  (5/dt)*(gbai*iau  +  gbbi*ibu  +  gbci*icu); 
vccc  =  -  (.5/dt)*(gcai*iau  +  gcbi*ibu  +  gcci*icu); 


vaf  =  vca  +  vcac; 
vbf  =  vcb  +  vcbc, 
vcf  =  vcc  +  vccc; 

err  =  (vaf-ova)*(vaf-ova)+(vbf-ovb)*(vbf-ovb)+(vcf-ovc)*(vcf-ovc); 

setva(vaf); 
set_vb(vbf), 
set_vc(vcf); 

//        cout  «  "Currents  "  «  iau  «  "  "  «  ibu  «  "  "  «  icu  «  "\n"; 
//        cout  «  "Voltages  "  «  vca  «  "  "  «  vcb  «  "  "  «  vcc  «  "\n"; 
//        cout  «  "I  Corr  V  "  «  vcac  «  "  "  «  vcbc  «  "  "  «  vccc  «  "\n"; 
//        cout  «  "Error  =  "  «  err  «  "\n", 
return(err), 
} 

void  bus:  :report() 
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{ 

cout  «  "Bus  "  «  nodeno  «  "  Has  "  «  nlines  «  "  Lines\n", 
for  (int  j=0;  j<nlines;  j++) 
cout  «  "Line  "  « j  « 
"  Between  Nodes  "  «  lptr[j]->report_node_a()  «  "  And  "  « 
lptr[j]->report_node_b()  «  "  Line  Pointer  "  «  &lptr[j]  «  "\n", 
} 

void  vbus::report() 

{ 
cout  «  "Voltage  Bus  V=  "  «  vamp  «  "  Phase=  "  «  vphase  «  "\n"; 

} 

void  line::setline  (int  na,  int  nb,  double  x,  double  xz,  double  ra, 
double  omza=377, 
int  saa  =  1 ,  int  sba  =  1 ,  int  sea  =  1 , 
double  iaa  =  0,  double  iba  =  0,  double  ica  =  0) 

{ 

nodea  =  na, 
nodeb  =  nb; 
xs  =  (2.0*x  +  xz)/3.0, 
xm  =  (xz  -  x)/3.0; 
omz  =  omza; 
r  =  ra; 
sa  =  saa; 
sb  =  sba, 
sc  =  sea; 
iao  =  ia  =  iaa; 
ibo  =  ib  =  iba; 
ico  =  ic  =  ica; 

cout  «  "Line:  setline  xs  =  "  «  xs  «  "  xm  =  "  «  xm  «  "  r  =  "  «  r 
«  "\n  Sw  =  "  «  sa  «  sb  «  sc  «  "  i  =  "  «  ia  «  "  "  «  ib 
«  "  "  «  ic  «  "\n\n"; 

} 

void  line::set_bus_pointers(bbus*  busa,  bbus*  busb) 

{ 

cout  «  "line::set_bus_pointers  "  «  busa  «  "  "  «  busb  «  "\n", 

busa  =  busa, 

busb  =  busb, 
} 

void  line::init_currents  (double  iaa,  double  iba,  double  ica) 

{ 
iao  =  ia  =  iaa; 
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ibo  =  ib  =  iba, 
ico  =  ic  =  ica; 

cout  «  "Line  Init  Currents  "  «  ia  «  "  "  «  ib  «  "  "  «  ic  «  "\n", 
} 

void  line::setup  ()  //  Conductance  Parameters 

{ 

double  d,  gs,  gm; 
//  int  abs(int); 

if  (abs(sa)  +  abs(sb)  +  abs(sc)  ==  3)    //  all  in 

{ 

d  =  xs  *  xs  -  2.0  *  xm  *  xm  +  xs  *  xm, 

gs  =  omz  *  (xs  +  xm)  /d; 

gm  =  -  omz  *  xm  /  d, 

gaa  =  gbb  =  gcc  =  gs; 

gab  =  gba  =  gac  =  gca  =  gcb  =  gbc  =  gm; 

} 
else  if  (abs  (sa)  +  abs  (sb)  +  abs  (sc)  ==  2)  //  one  line  out 

{ 

d  =  xs  *  xs  -  xm  *  xm, 

gs  =  omz  *  xs  /  d; 

gm  =  -omz  *  xm  /  d; 

if  (sa  ==  0)  //  line  A  out 

{ 

gaa  =  gab  =  gac  =  gba  =  gca  =  0, 

gbb  =  gcc  =  gs; 

gbc  =  gcb  =  gm, 
} 

else  if  (sb  ==  0)  //  line  B  out 

{ 

gbb  =  gab  =  gba  =  gbc  =  gcb  =  0, 
gaa  =  gcc  =  gs, 
gac  =  gca  =  gm; 

} 

else  if  (sc  =  0) 

{ 

gcc  =  gca  =  gac  =  gbc  =  gcb  =  0, 
gaa  =  gbb  =  gs; 
gab  =  gba  =  gm, 

} 
else  {  cout  «  "Blew  One  Line  Out  Case!\\n",} 

} 
else  if((abs(sa)  +  abs(sb)  +  abs(sc))  ==  1) 
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if(abs(sa)  =  l) 

{ 

gaa  =  omz/xs; 
gab  =  gba  =  gbb  =  gac  =  gca  =  gcb  =  gbc  =  gcc  =  0; 

} 
else  if  (abs(sb)  ==  1) 

I 

gbb  =  omz/xs; 
gaa  =  gab  =  gba  =  gbc  =  gcb  =  gca  =  gac  =  gcc  =  0, 

} 
else  if  (abs(sc)  ==  1) 

{ 

gcc  =  omz/xs; 
gaa  =  gab  =  gba  =  gbb  =  gac  =  gca  =  gbc  =  gcb  =  0; 

} 
else 


cout  «  "Blew  One  Line  IN  case!  \n"; 


} 
else  if((abs(sa)  +  abs(sb)  +  abs(sc))  ==  0) 

{ 

gaa  =  gab  =  gba  =  gca  =  gac  =  gbb  =  gbc  =  gcb  =  gcc  =  0; 

} 
else  {cout  «  "Bad  Combination  of  Switch  States!\n",} 


i 


cout  «  "Line  :Setup()  G  =  "  «  gaa  «  "  "  «  gab  «  "  "  «  gac 

«  "\n  "  «  gba  «  "  "  «  gbb  «  "  "  «  gbc 

«  "\n  "  «  gca  «  "  "  «  gcb  «  "  "  «  gcc  «  "\n\n"; 


void  line 
void  line 
void  line 


open_a()  { sa  = 
open_b()  { sb  = 
opencQ  {sc  = 


i;} 
■i;} 
■i;} 


void  line 
void  line 
void  line 


close_a()  {sa=  1;} 
close_b()  {sb  =  1;} 
close_c()  {sc=  1;} 


int  linexurrentmonitor  ()  //  Checks  for  zero  crossings 

I 
int  ip  =  0; 

if  ((sa  ==  -1)  &&  (iao  !=0)  &&  (iao  *  ia  <0))  ip  +=  1; 
if  ((sb  =  -1)  &&  (ibo  1=0)  &&  (ibo  *  ib  <0))  ip  +=  2; 
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if  ((sc  ==  -1)  &&  (ico  !=0)  &&  (ico  *  ic  <0))  ip  +=  4, 

//re-set  old  currents 

iao  =  ia, 
ibo  =  ib; 

ico  =  ic; 

return  (ip), 
} 

void  line::re_setup  (int  ip)    //  re-build  line  model 

{ 

if  (ip  &  1)  {  sa  =  0;  ia  =  0;}  //  Phase  A  opening 
if  (ip  &  2)  {  sb  =  0;  ib  =  0,}  //  Phase  B  opening 
if  (ip  &  4)  {  sc  =  0;  ic  =  0,}  //  Phase  C  opening 
setup(); 

} 

void  line::rkl  (double  dt) 

{ 

double  va,  vb,  vc; 

va  =  busa->a_voltage()  -  busb->a_voltage(), 
vb  =  busa->b_voltage()  -  busb->b_voltage(), 
vc  =  busa->c_voltage()  -  busb->c_voltage(); 

kla  =  gaa  *  (va  -  r  *  ia)  *  dt 
+  gab  *  (vb  -  r  *  ib)  *  dt 
+  gac  *  (vc  -  r  *  ic)  *  dt; 

klb  =  gba  *  (va  -  r  *  ia)  *  dt 
+  gbb  *  (vb  -  r  *  ib)  *  dt 
+  gbc  *  (vc  -  r  *  ic)  *  dt; 

klc  =  gca  *  (va  -  r  *  ia)  *  dt 
+  gcb  *  (vb  -  r  *  ib)  *  dt 
+  gcc  *  (vc  -  r  *  ic)  *  dt; 

iap  =  ia  +  5*kla,  //usable  for  next  rk  step,  both  internal  and 

ibp  =  ib  +  5*klb;  //outside,  reported  by  i_xp() 

icp  =  ic  +  5*klc, 
} 

void  line::rk2  (double  dt) 
{ 
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double  va,  vb,  vc; 

va  =  busa->a_voltage()  -  busb->a_voltage(); 
vb  =  busa->b_voltage()  -  busb->b_voltage(), 
vc  =  busa->c_voltage()  -  busb->c_voltage(); 

k2a  =  gaa  *  (va  -  r  *  iap)  *  dt 
+  gab  *  (vb  -  r  *  ibp)  *  dt 


+  gac  * 

(vc  - 

r  *  icp)  * 

dt; 

k2b  =  gba 

*  (va 

- r  *  iap) 

*dt 

+  gbb* 

(vb- 

r  *  ibp)  * 

dt 

+  gbc  * 

(vc  - 

r  *  icp)  * 

dt, 

k2c  =  gca 

*  (va 

- r  *  iap) 

*dt 

+  gcb  * 

(vb- 

r  *  ibp)  * 

dt 

+  gcc  * 

(vc  - 

r  *  icp)  * 

dt; 

iap  =  ia  +  5*k2a;  //usable  for  next  rk  step,  both  internal  and 
ibp  =  ib  +  .5*k2b;  //outside,  reported  by  i_xp() 
icp  =  ic  +  5*k2c; 
} 

void  line::rk3  (double  dt) 

{ 

double  va,  vb,  vc; 

va  =  busa->a_voltage()  -  busb->a_voltage(), 
vb  =  busa->b_voltage()  -  busb->b_voltage(); 
vc  =  busa->c_voltage()  -  busb->c_voltage(), 

k3a  =  gaa  *  (va  -  r  *  iap)  *  dt 
+  gab  *  (vb  -  r  *  ibp)  *  dt 
+  gac  *  (vc  -  r  *  icp)  *  dt; 

k3b  =  gba  *  (va  -  r  *  iap)  *  dt 
+  gbb  *  (vb  -  r  *  ibp)  *  dt 
+  gbc  *  (vc  -  r  *  icp)  *  dt; 

k3c  =  gca  *  (va  -  r  *  iap)  *  dt 
+  gcb  *  (vb  -  r  *  ibp)  *  dt 
+  gcc  *  (vc  -  r  *  icp)  *  dt; 

iap  =  ia  +  k3a;  //usable  for  next  rk  step,  both  internal  and 
ibp  =  ib  +  k3b,  //outside,  reported  by  i_xp() 
icp  =  ic  +  k3c; 

} 
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void  line::rk4  (double  dt) 

{ 

double  va,  vb,  vc; 

va  =  busa->a_voltage()  -  busb->a_voltage(); 
vb  =  busa->b_voltage()  -  busb->b_voltage(); 
vc  =  busa->c_voltage()  -  busb->c_voltage(); 

k4a  =  gaa  *  (va  -  r  *  iap)  *  dt 
+  gab  *  (vb  -  r  *  ibp)  *  dt 
+  gac  *  (vc  -  r  *  icp)  *  dt; 

k4b  =  gba  *  (va  -  r  *  iap)  *  dt 
+  gbb  *  (vb  -  r  *  ibp)  *  dt 
+  gbc  *  (vc  -  r  *  icp)  *  dt; 

k4c  =  gca  *  (va  -  r  *  iap)  *  dt 
+  gcb  *  (vb  -  r  *  ibp)  *  dt 
+  gcc  *  (vc  -  r  *  icp)  *  dt; 
} 

void  line::rk  ()  //  Finishes  the  Runge-Kutta 

{ 

ia  +=  (kla  +  2.0  *  k2a  +  2.0  *  k3a  +  k4a)  /  6.0, 

ib  +=  (klb  +  2.0  *  k2b  +  2.0  *  k3b  +  k4b)  /  6.0; 

ic  +=  (klc  +  2.0  *  k2c  +  2.0  *  k3c  +  k4c)  /  6.0; 

} 

void  line::report() 

{ 

cout  «  "Line  From  "  «  nodea  «  "  To  "  «  nodeb  « 
..  p=  -  «  r  «  ■■  xs  =  ••  «  xs  «  ••  xm  =  »  «  xm  «  "\n»; 

cout  «  "Bus  Pointers  are  "  «  busa  «  "  "  «  busb  «  "\n", 

} 

void  line::report(char*  label) 

{ 

cout  «  "line:  report  "  «  label  «  "\n", 

cout  «  "Line  From  "  «  nodea  «  "  To  "  «  nodeb  « 
"  r=  "  «  r  «  "  xs  =  "  «  xs  «  "  xm  =  "  «  xm  «  "\n"; 

cout  «  "Bus  Pointers  are  "  «  busa  «  "  "  «  busb  «  "\n", 
} 

int  line  :report_node_a()  {return(nodea),} 
int  line::report_node_b()  {return(nodeb),} 


135 


double  line::va_other_end  (int  thisnode) 

{ 
//     cout  «  "line.:va_other_end  (  '*«  thisnode  «  ")\nM; 
//     cout  «  "nodea  =  "  «  nodea  «  "  nodeb  =  "  «  nodeb  «  "\n", 
//     cout  «  "busa  pointer  =  "  «  busa  «  "  busb  pointer  "  «  busb  «  "\n"; 
//     cout  «  "busa  v  =  "  «  busa->a_voltage()  « 
//       "  busb  v  =  "  «  busb->a_voltage()  «  "\n"; 
if  (thisnode  ==  nodea) 

return(busb->a_voltage()); 
else  return(busa->a_voltage()); 
} 

double  line::vb_other_end  (int  thisnode) 

{ 

if  (thisnode  ==  nodea) 

return(bu  sb->b_voltage( )) ; 
else  return(busa->b_voltage()); 

} 

double  line::vc_other_end  (int  thisnode) 

{ 

if  (thisnode  ==  nodea) 
return(busb->c_voltage()), 

else  return(busa->c_voltage()), 
} 

int  line::report_other_node  (int  thisnode) 

{ 

if  (thisnode  ==  nodea) 

return(nodeb), 
else  return(nodea); 
} 

double  line  :ias(int  thisnode) 

{ 

if  (thisnode  =  nodea) 
return(ia); 

else  return(-ia), 
} 

double  line:  :ibs(int  thisnode) 

{ 

if  (thisnode  ==  nodea) 
return(ib), 
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else  return(-ib), 
} 

double  line::ics(int  thisnode) 

{ 

if  (thisnode  ==  nodea) 
return(ic); 

else  return(-ic); 
} 
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